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Executive  Summary 

The  objective  of  this  report  is  to: 

(1)  present  an  algorithm  for  regridding  and  transformation  of  dig¬ 
ital  terrain  elevation  matrices  that  will  efficiently  run  on  a 
minicomputer , 

(2)  benchmark  transformation  and  I/O  times  on  four  candidate  mini¬ 
computer  systems, 

(3)  make  recommendations  as  to  which  minicomputer  systems  will  per¬ 
form  most  efficiently,  and 

(4)  recommend  modifications  to  DMAHTC  current  Mosaicking  and  Re¬ 
gridding  Software  System  (MARS) . 

The  motivation  for  the  study  effort  came  from  the  realization  that 
the  MARS  system,  which  is  implemented  on  a  UNIVAC  1100/81  mainframe,  has 
bottlenecks  that  result  from  certain  operating  system  constraints.  In 
looking  for  a  solution  to  this  problem  it  became  apparent  that  state-of- 
the-art  minicomputer  systems,  particularly  the  32  bit  machines,  could 
perform  as  well  or  better  than  the  UNIVAC  1100/81,  especially  if  the  re¬ 
gridding  and  transformation  functions  were  performed  in  a  sequential  manner 
rather  than  the  current  random  manner.  A  new  algorithm  was  designed  with  a 
view  toward  installing  it  on  a  "fast"  minicomputer.  The  new  algorithm  could 
be  implemented  on  the  UNIVAC  1100/81,  however,  it  would  not  achieve  the 
same  efficiency  as  it  will  on  a  suitable  minicomputer. 

Section  2  describes  the  new  algorithm  and  its  refinement.  The  new 


algorithm  rearranges  the  regridding  and  transformation  task  so  that  se¬ 
quential  accessing  methods  to  mass  storage  devices  can  replace  the  current 


random  access  methods.  Due  to  the  resulting  simplification  of  the  processing 
needs  the  new  algorithm  can  be  installed  on  a  minicomputer,  allowing  easy 
and  inexpensive  expansion  should  there  be  a  future  increase  in  workload. 

Section  3  examines  the  specifications  of  minicomputer  systems  which 
could  perform  the  new  algorithm  efficiently  enough  to  compete  with  the  UNIVAC 
1100/81.  Both  16  bit  and  32  bit  minicomputers  are  examined.  A  comparison  of 
relative  costs  is  made.  Special  hardware  options  which  can  significantly 
improve  execution  times  are  noted. 

Benchmark  tests  were  run  on  a  subset  of  the  minicomputers  presented  in 
Section  3  in  order  to  determine  if  preliminary  estimates  of  their  speeds  were 
correct.  In  Section  4  the  results  of  these  tests  are  given.  I/O  transmission 
speeds  are  shown  to  be  more  than  sufficient  for  the  implementation  of  the 
new  algorithm.  The  speed  at  which  a  minicomputer  could  perform  a  compute- 
bound  point-to-point  transformation,  necessary  for  the  regridding  task, 
shows  that  the  fastest  (and  most  expensive)  minicomputers  would  be  necessary 
if  each  point  in  the  Digital  Terrain  Matrix  were  to  be  rigorously  trans¬ 
formed.  However,  only  a  small  portion  of  the  points  are  rigorously  trans¬ 
formed  by  the  current  software  used  by  DMA.  Under  these  circumstances,  all 
the  minicomputers  examined  can  be  used  for  the  new  algorithm. 

Scetion  5  presents  an  overview  of  the  study,  conclusions,  and  recom¬ 
mendations.  General  conclusions  regarding  the  transportation  of  the  software 
from  the  UNIVAC  1100/81  to  a  minicomputer,  including  estimates  of  costs,  are 
given.  The  recommendation  is  made  that  a  32  bit  minicomputer  (PERKIN-ELMER 
3240)  be  chosen  for  the  implementation  of  the  new  algorithm.  This  recommen¬ 
dation  is  based  on  the  ease  of  implementation  possible  with  the  larger  address 
space  of  a  32  bit  computer  when  compared  to  that  of  a  16  bit  computer  and 
the  capability  of  a  32  bit  computer  in  handling  future  DMA  needs. 


Table  of  Contents 


PAGE 

Executive  Summary . iv 

1.  Introduction . 1-1 

2.  Algorithms . 2-1 

2.1  The  Proposed  Algorithm . 2-1 

2.2  The  Refined  Algorithm . 2-3 

2.2.1  Notation . 2-3 

2.2.2  Regridding  and  Binning . 2-3 

2.2.3  Interpolation . 2-6 

2.2.4  Coordinate  Coding . 2-6 

2.2.5  Sampling  Grid  Generation . 2-6 

2.2.6  Implementation . 2-8 

2.3  Further  Analysis  of  the  Refined  Algorithm . 2-8 

2.3.1  Assumptions  Concerning  Hardware . 2-9 

3.  Hardware . 3-1 

3.1  The  DATA  GENERAL  Eclipse  S/250 . 3-3 

3.2  The  PERKIN-ELMER  3220  and  3240 . 3-5 

3.3  Comparison  of  Systems . 3-8 

4.  Benchmarks . 4-1 

5.  Conclusions  and  Recommendations . 5-1 

5.1  General  Conclusions . 5-1 

5.2  Conclusions  Regarding  Coordinate  Transformations. ...  5-2 

5.3  Final  Conclusions . 5-3 


Appendix  A  —  Study  of  Transformation  Package 

Appendix  B  —  Formal  Minicomputer  Manufacturer's  Quotations 


vi 


List  of  Figures 


Figure  1.  Flow  diagram  for  the  initially  proposed 
algorithm 

Figure  2.  Sample  Input  and  Output  of  Initial 
MiniMARS  Test 

Figure  3.  All  output  grid  points  will  lie  in  one 
of  N  +  1  aisles  for  an  NxN  input  DTM 

Figure  4.  Elevations  of  transformed  output  points 
are  interpolated  from  elevations  of 
input  points 

Figure  5.  Flow  diagram  for  the  refined  MiniMARS 
algorithm 

Figure  6.  Peripheral  flow  diagram  for  the  refined 
MiniMARS  algorithm 


vii 


List  of  Tables 


PAGE 

Table  1.  Comparative  Prices . 3-7 

Table  2.  Transformation  Benchmarks  for  MiniMARS . 4-2 

Table  3.  Geographic  to  UTM  Times  on  Four  Minicomputers ....  4-2 

Table  4.  Random  Reads  from  Disk . 4-4 

Table  5.  Sequential  Reads  from  Disk . 4-4 


viii 


INTRODUCTION 

The  Defense  Mapping  Agency  Hydrographic/Topographic  Center  (DMAHTC)  has 
a  requirement  to  convert  a  large  photogrammetr ic  source  database  to  digital 
terrain  elevation  matrices.  This  conversion  is  currently  achieved  in  large 
part  by  scanning  raw  photogrammetric  data  with  specialized  hardware  (Univer¬ 
sal  Automatic  Map  Compilation  Equipment  or  UNAMACE) ,  the  output  of  which  is 

a  matrix  of  elevation  values  (Digital  Terrain  Matrix  or  DTM) .  DTM's  are 

0 

combined  to  form  1  squares  and  stored  in  a  digital  database  using  DMAHTC 
Standard  Format.  The  process  of  combining,  called  mosaicking,  involves  the 
regridding  and  transformation  of  the  DTM's.  This  is  currently  accomplished 
using  a  complex  system  of  software  called  the  Mosaicking  And  Regridding 
System  (MARS)  developed  and  implemented  at  DMAH1C  by  L.N.K.  Corporation  in 
1978. 

MARS,  designed  to  run  in  batch  mode  on  a  UNIVAC  1100/81  computer,  has 
been  vital  to  DMA  in  handling  the  conversion  of  UNAMACE  data.  However, 
after  working  at  DMAHTC  with  the  system  on  a  daily  basis  for  two  years  and 
gaining  an  intimate  knowledge  of  the  operator  interface,  machine  resource 
requirements  and  wall  clock  time  necessary  for  the  regridding  and  transfor¬ 
mation  process,  L.N.K.  determined  that  some  portions  of  the  MARS  system  were 
in  effect  bottlenecks.  These  bottlenecks  are  due  to  constraints  imposed  by 
the  UNIVAC  operating  system.  In  particular,  the  limiting  of  the  memory  address 
space  leads  to  a  fragmentation  of  the  data  being  processed,  thus  excessive 
random  accesses  to  mass  storage  devices  are  required. 

After  collection  of  statistics  on  the  regridding  and  transformation 
process  and  upon  closer  examination  of  both  the  details  of  the  process  and 
the  capabilities  of  state-of-the-art  minicomputers,  L.N.K.  estimated  that 
the  regridding  and  transformation  task  could  be  performed  on  a  dedicated 
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minicomputer  more  efficiently  than  on  the  UNIVAC  1100/81.  In  addition,  the 
implementation  of  a  "MiniMARS”  would  allow  for  easy  expansion  to  handle  any 
increase  in  processing  needs  simply  by  adding  to  the  number  of  minicomputers 
used  for  the  task. 

L.N.K.  therefore  proposed  a  study  effort  to: 

(1)  Develop  a  new  algorithm  for  the  regridding  and  transformation  pro¬ 
cess  which  would  significantly  decrease  the  number  of  random  accesses 
to  mass  storage  devices, 

(2)  Examine  a  selection  of  minicomputer  systems  and  prove  they  could 
efficiently  and  effectively  perform  the  new  algorithm, 

(3)  Recommend  DMAHTC  purchase  any  one  of  four  select  systems,  and 

(4)  Recommend  modifications  to  the  current  MARS  system  which  could 
improve  throughput  on  the  UNIVAC  1100/81.  These  modifications  are 
not  an  implementation  of  the  new  algorithm,  but  an  optimization  of 
the  code  used  in  the  current  MARS. 

In  the  study  benchmark  tests  were  devised  and  implemented  on  a  number 
of  minicomputer  systems  to  determine  possible  transformation  and  1/0  times 
for  the  regridding  and  transformation  task.  Using  the  results  from  these  tests 
L.N.K.  confirmed  that  state-of-the-art  minicomputers  (PERKIN-ELMER  3240,  DATA 
GENERAL  ECLIPSE  S/250,  PDP  11/44)  can  compete  with  large  mainframes  such  as 
the  UNIVAC  1100/81  when  the  minicomputers  are  dedicated  to  the  task  at  hand. 
The  details  of  the  algorithm,  tests,  and  recommendations  are  included  in  this 
report. 


1-2 


2. 


ALGORITHMS 


2 . I  The  Proposed  Algorithm 

Universally,  the  regridding  and  transformation  task  consists  of 

(1)  Definition  and  creation  of  a  target  coordinate  system  grid, 

(2)  Input  and  transformation  of  input  digital  elevation  data,  and 

(3)  Creation  of  the  output  model. 

Regridding  and  transformation,  normally  performed  using  random  access  methods 
with  the  MARS  software  system,  has  been  restructured  in  the  new  algorithm 
into  a  series  of  sequential  file  processing  steps.  Each  step,  due  to  simpli¬ 
city,  can  be  efficiently  implemented  on  a  minicomputer  system. 

The  proposed  new  algorithm  is  a  4  step  process.  Initially  the  target 
coordinate  system  grid  is  defined.  Each  point  in  this  grid  (x',y')  is  trans¬ 
formed  to  the  input  model's  coordinate  system,  thus  producing  a  set  of 
quadruples  (x’,y',x,y).  These  quadruples  are  then  sorted  according  to  the 
(x,y)  coordinates.  The  sorted  quadruples  are  now  merged  with  the  input  model, 
interpolating  elevation  values  for  (x,y)  as  is  currently  done  with  the  MARS 
system.  The  result  of  this  step  is  a  set  of  triples  (x',y',z')  which  are  in 
the  input  model's  order.  Therefore  these  triples  are  resorted  according  to 
the  (x',y')  coordinates,  producing  a  file  with  the  correct  order  for  the 
output  model.  Figure  1  illustrates  the  steps  of  this  process.  Note  that  the 
transformation  from  one  coordinate  system  to  the  other  is  considered  to  be 
done  by  a  "black  box",  i.e.  any  transformation  function  can  be  used  with 
this  system,  including  those  used  by  MARS. 

An  elementary  version  of  the  algorithm  was  implemented  on  a  UNIVAC  1108 
computer.  Initially  a  small  input  model  (20  x  20  points)  was  hand-created.  A 
target  coordinate  system  was  defined  as  a  small  rotation  and  translation  of  the 
input  coordinate  system.  A  file  of  400  quadruples  (x',y',x,y)  was  created. 
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Figure  1.  Flow  diagram  for  initially  proposed  algorithm 


The  quadruple  file  was  sorted  according  to  the  (x,y)  coordinates  using  the 
UNIVAC  system  sort  routine  SORTSDF.  The  file  of  quadruples  and  the  input 
model  were  then  merged,  interpolating  new  elevation  values  by  using  a  simple 
four  neighbor  average.  The  output  of  this  step  (x'.y'.z')  was  sorted  according 
to  the  (x',y')  coordinates  and  then  plotted.  Figure  2  shows  both  the  input 
model  and  the  resulting  output  model.  Numerical  elevation  values  are  repre¬ 
sented  by  alphanumeric  characters.  It  should  be  noted  that  in  a  real  system 
no  elevation  data  would  be  thrown  away.  In  this  case  data  was  clipped  only 
for  uniformity  of  plot  size. 

The  algorithm  thus  has  been  graphically  shown  to  be  correct.  However 
on  further  analysis  it  was  found  that  by  increasing  the  size  of  the  memory 
space  and  using  in-core  addressing  methods  the  sort  steps  could  be  dropped, 
greatly  increasing  the  throughput  of  the  system.  The  refined  algorithm  is 
detailed  below. 

2,2  The  Refined  Algorithm 

2.2.1  NOTATION 

As  in  the  proposed  system  primes  are  used  to  indicate  points  in  the  output 
grid  coordinates,  such  as  (x',y').  Unprimed  quantities  are  reserved  for  input 
grid  coordinates,  such  as  (x,y).  T(x',y')  or  T(x',y',z)  is  the  output  grid 
point  (x',y')  transformed  into  the  input  grid  coordinate  system. 

2.2.2  REGRIDDING  AND  BINNING 

As  Figure  3  shows,  output  grid  points  will  be  transformed  into  the 
input  Digital  Terrain  Matrix  (DTM)  coordinate  system  so  that  elevation  values 
can  be  determined  for  each  one  by  interpolation.  Assuming  an  N  x  N  input  DTM, 
every  transformed  grid  point  must  lie  in  one  of  N  +  1  "aisles"  or  bins.  For  a 
1024  x  1024  innut  DTM  there  are  only  1025  bins. 
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/  MODEL  AFTER  TRANSFORMATION 

/ 

//  TO  OUTPUT  COORDINATE  SYSTEM 

The  line  segments  at  right  angles 
enclose  the  same  region  in  the  two 
models.  Note  that  due  to  inter¬ 
polation  the  elevation  values  are 
slightly  changed  at  region  borders. 


Figure  2.  Sample  Input  and  Output  of  Initial  MiniMARS  Test 
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2.2.3  INTERPOLATION 

During  interpolation,  the  elevation  of  each  output  grid  point 
(0^,02* ...  .^28  FFSure  3)  will  be  determined  from  the  elevations  of  only 
4  points  in  the  input  DTM,  which  lie  on  only  two  profiles.  Figure  4  shows 
the  4  "neighbors"  used  to  compute  the  elevation  at  point  T(x',y').  If  the 
input  grid  points  are  never  more  than  128  meters  apart  at  true  scale  a 
fraction  of  6-bit  precision  for  T&',y')  will  yield  2  meter  positional 
accuracy  in  the  (x,y)  space. 

2.2.4  COORDINATE  CODING 

Using  16  bits  overall  and  6  bits  for  a  fraction,  T(x',y')  =  (x,y)  will 
be  coded  as  follows 


X=  Dx1Dx2---Dx10'Fx1Fx2'--Fx6 

y~  DylDy2-*-Dy10’FylFy2-.-Fy6 

The  consequences  of  2  meter  accuracy  at  up  to  128  meter  spacing  is  a  limit 
on  the  model  size  to  1024  x  1024  if  the  above  16-bit  coordinates  are  to  be 
used.  However,  the  above  bit  pattern  can  be  varied  as  is  needed.  Using 
fewer  fraction  bits  allows  a  larger  DTM  to  be  manipulated. 


2.2.5  SAMPLING  GRID  GENERATION 

The  output  grid  points  will  be  generated  and  transformed  into  the  j 
input  coordinate  space.  During  generation  each  point  will  be  represented  | 
by  a  quadruple  (x’,y',x,y),  where  (x,y)  -  T(x',y'),  which  will  immediately  | 
be  placed  in  the  correct  aisle  for  later  sequential  merging  with  the  input  I 


2-6 


ELEVATION  AT  T(x',y')  IS  INTERPOLATED  FROM  THE 
ELEVATIONS  AT  (X^Yj.),  (X^Y^,  (X1,Y2),  AND 
(X2,Y2).  IF  INPUT  SPACING  IS  LESS  THAN  OR  EQUAL 
TO  128  METERS,  A  6  BIT  FRACTION  WILL  YIELD  2 
METER  ACCURACY. 


Figure  4.  Elevations  of  transformed  output  points  are 

interpolated  from  elevations  of  input  points. 


DTM.  The  rough  algorithm  is  as  follows: 

(Al)  Get  input  DTM  specs  from  the  header  and  define  the  output  grid 
and  transformation  T. 

(A2)  Generate  and  bin  each  point  of  the  output  grid 

(A2.1)  For  each  point  (x',y')  produce  the  quadruple  (x',y',x,y) 

where  (x,y)  =  T(x',y')  and  where  x  ^  D  ]D  1...D  10'F  IF  2-.-F 

X  X  x 

and  y  =  DylDy2...Dyl0-FylFy2-..Fy6 
(A2.2)  Place  each  quadruple  (x',y',x,y)  into  aisle  A  which  is 
DxlDx2---  Dx10  in  binary. 

2.2.6  IMPLEMENTATION 

No  aisle  can  receive  more  than  1024  / 2  quadruples  (the  square  root 

value  comes  from  the  diagonal  of  the  matrix)  hence  8  times  this  (about 

15000)  bytes  of  storage  are  needed  for  each  aisle.  In  main  memory  we  can 

statically  allocate  1025  blockr  of  size  25  quadruples  (i.e.  200  bytes)  each. 

When  full,  each  25  quadruple  block  will  be  output  to  a  known  track  storing 

20 

the  appropriate  aisle.  There  will  thus  be  about  2  /25  writes  (about  40000). 

23  5 

This  will  require  40000  x  43  ms  average  access  time  plus  2  bytes/ (8x10 
bytes/second  transfer  rate)  or  about  1730  seconds,  approximately  1/2  hour. 

A  refined  algorithm  will  now  be  given  which  will  reduce  the  number  of  disk 
writes,  and  hence  the  execution  time,  by  half. 

2 . 3  Further  Analysis  of  the  Refined  Algorithm 

An  in  depth  analysis  of  the  refined  algorithm  using  data  on  off-the-shelf 
minicomputer  systems  is  given  on  the  following  pages.  Though  the  DATA  GENERAL 
S/250  was  used  for  the  basis  of  the  times  given,  all  the  machines  examined 
in  Section  3  have  similar  1/0  speeds.  Therefore  the  use  of  the  DATA  GENERAL 
times  should  not  be  construed  as  a  recommendation  of  the  S/250. 
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2.3.1  ASSUMPTIONS  CONCERNING  HARDWARE 

The  refined  algorithm  which  follows  refers  to  general  specifications 

for  high  speed  minicomputer  systems,  for  example  the  DATA  GENERAL  S/250  or 

the  PERKIN-ELMER  3240.  Memory  size  necessary  for  efficient  implementation 

1 8 

of  the  algorithm  should  be  at  least  2  bytes  (256K  bytes) .  Peripheral 
equipment  also  needed  include  two  independent  disk  drives  (i.e.,  capable 
of  performing  operations  simultaneously)  and  two  1600BPI  9-track  tape  drives. 
The  disk  drives  assumed  each  have  five  surfaces  and  815  cylinders.  Each  track 
(intersection  of  one  surface  with  a  cylinder)  contains  24  sectors  of  512 
bytes.  The  average  access  time  of  the  disk  drives  is  an  important  parameter 
and  should  be  close  to  the  43  millisecond  access  time  of  the  DATA  GENERAL 
drives.  The  tape  drives,  though  not  as  important  with  respect  to  speed,  are 
assumed  to  be  capable  of  reading  and  writing  at  75  inches/second. 
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Step  1.  Generate  output  grid,  transform,  and  "aisje  sort"  to  disk. 
Algorithm 

Read  header  of  input  DTM  and  define  output  grid. 

Generate  and  transform  output  grid  points  (x',y')  to  get  quads 
(x',y',x,y)  with  (x,y)  =  T(x',y'). 

Output  quads  in  blocks  of  512  bytes,  or  64  quads,  to  one  of  256 
cylinders  on  disk  according  to  first  8  bits  of  x.  (i.e.  each 
cylinder  would  store  4  aisles  of  data) 

Main  Memory  Data  Structures 

256  block  buffers  @  512  bytes  each  (128K  bytes) 

3  output  buffers  @  512  bytes  each  (1.5K  bytes) 

256  block  headers  @  4  bytes  each  (ik  bytes) 

Disk  Structure 

bin  K  stored  on  cylinder  K,  0  _<  K  _<  255 
K  selected  as  first  8  bits  of  x  in  quad 

60,000  bytes  per  cylinder  allows  7,000  quads  max  per  4  aisles 

n  0 

Timing  for  2  points  assuming  no  1/0  overlap 

Transform:  3  to  45  minutes 
Output:  8  Mbvtes/512  =  16K  blocks 

16K  blocks  x  50  x  10  seconds  each  =  800  seconds  =  13  minutes 

Step  2.  Collate  quads  with  input  DTM  profiles  and  sort  on  x'  on  output 
Algorithm 

Read  first  5  input  DTM  profiles. 

Read  next  block  of  quads  in  one  of  4  aisles  between  the  5  profiles. 


2-10 


For  each  quad  (x',y',x,y)  use  (x,y)  to  interpolate  z'.  Then  place 
(x',y',z')  in  output  block  buffer  according  to  the  first  byte  of  x'. 
When  the  current  4  aisles  are  exhausted  move  on  to  the  next  4  input 
DTM  profiles. 

Main  Memory  Data  Structures 

Tape:  10  input  profile  buffers  @  2,050  bytes  each  (20K  bytes) 

256  block  buffers  @  512  bytes  each  (128K  bytes) 

256  block  headers  @  4  bytes  each  (IK  bytes) 

Disk  out:  3  output  buffers  @  512  bytes  each  (1.5K  bytes) 

Disk  in:  3  input  buffers  @  2,048  bytes  each  (6K  bytes) 

Disk  Structure 

Input:  quad  file  as  in  Step  1  output 

Output  triples 

bin  K  stored  on  cylinder  K,  0  _<  K  <_  255 
K  selected  as  first  byte  of  x' 

85  triples  per  512  byte  block,  hence  10,200  triples  per  cylinder 
this  allows  10  profiles  @  1,020  points  each;  typically  only 
4  will  exist  since  x'  has  only  10  bits  of  integer 

on 

Timing  for  2  _  points  assuming  no  I/O  overlap 
input  from  tape  1000  blocks  @  2"  each 
2000"  @  75ips  =  30  seconds 
input  from  diskj 

4K  blocks  (of  2048  bytes)  x  50  x  10”^  seconds  each  =  200  seconds 
output  to  disk2 

12K  blocks  (of  512  bytes)  x  50  x  10”^  seconds  =  600  seconds 
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main  memory  processing  (estimated  interpolation  time) 


=  180  seconds 

TOTAL  TIME  ~  1,000  seconds 

Final  Model  Output 
Algorithm 

Read  entire  next  bin  (cylinder)  of  triples  and  assemble  up  to  10 
profiles  of  z'  values  by  selecting  the  profile  with  x'  and  the  point 
position  with  y' . 

Repeat  for  all  bins. 


t 

X 


Place  z'  in  position  addressed  by  x'  and  y' 

Ma in  Memory  Data  Structures 
output  buffers  for  output  tape 

20  profiles  x  1,500  points  x  2  bytes  each  =  60,000  bvtes 
for  2  sets  of  buffers 
input  buffers  tor  triples 

2  buffers  0  2,048  bytes  each 


r  i 


20 

Timing  assuming  no  overlap  for  2  points 
writing  output  tape  =  30  seconds 

in  core  =  200  seconds  (very  conservative  estimated  time) 
input  of  3K  blocks  of  2,048  bytes  each 

3K  x  50  x  10  ^  seconds  =  150  seconds 
TOTAL  TIME  =  380  seconds 

TOTAL  TIMES  FOR  THE  3  STEPS 


Step  1. 

17  to  58  minutes* 

(clear  bottleneck) 

Step  2. 

16  minutes 

Step  3. 

7  minutes 

TOTAL 

40  to  80  minutes 

*high  value  for  rigorous  point-to-point  transformation;  low  value  for  maximum 
interpolation 

Figures  5  and  6  illustrate  the  refined  algorithm.  The  sort  steps  have 

essentially  been  eliminated  through  the  use  of  direct  addressing  methods 

(both  ..lain  memory  and  disk).  The  refined  algorithm  requires  a  larger  (256K  = 
18 

2  )  memory  space  than  is  normally  directly  addressable  by  a  16  bit  mini¬ 

computer. 
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Figure  5.  Flow  diagram  for  the  refined  MiniMARS  algorithm 
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Figure  (.  Peripheral  flow  diagram  for  the  refined  MiniMARS  algorithm 
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HARDWARE 


The  refined  regridding  and  mosaicking  process  requires  a  fair  amount 
of  main  memory  (256K  bytes)  in  order  for  it  to  perform  fast  enough  to  be 
more  efficient  than  the  methods  used  in  the  MARS  software  system  on  the 
UNIVAC  1108.  The  basic  requirements  of  the  minicomputer  hardware  system 
are: 

1)  256K  byte  main  memory  (minimum) 

2)  Floating  point  instructions  (hardware,  not  software  interpreted) 

3)  2  fast  disk  drives  (average  a~cess  time  approximately  45 

milliseconds),  capacity  of  each  drive  at  least  50M  bytes 

4)  Two  1600  bpi  9  track  tape  drives,  75  inch  per  second  read/write 

speed  (or  faster) . 

In  addition,  some  method  of  producing  hardcopy  output  for  monitoring 
the  regridding  and  mosaicking  process,  such  as  a  printing  console  terminal, 
and  software  comprising  a  FORTRAN  compiler  and  text  editor  are  essential. 

Table  1  gives  a  price  breakdown  of  4  off-the-shelf  minicomputer  systems 
which  meet  these  requirements.  The  prices  stated  are  based  on  the  formal 
manufacturer's  quotes  contained  in  Appendix  B.  The  formal  quotes,  compiled 
by  sales  representatives  of  the  companies  DATA  GENERAL,  PERKIN-ELMER,  and 
DIGITAL  EQUIPMENT  contain  certain  discrepencies  in  system  configuration 
which  are  outlined  on  the  quotation  sheets.  The  discrepencies  can  increase 
the  bottom  line  costs  for  the  systems  by  as  much  as  20%.  Table  1  should  be 
used  as  the  cost  guideline  in  place  of  the  Quotation  sheets.  All  prices  in 
Table  1  reflect  the  manufacturer's  GSA  discount,  15%  for  DATA  GENERAL  and 
DIGITAL  EQUIPMENT,  21%  for  PERKIN-ELMER. 

The  systems  examined  were  felt  to  be  adequate  for  the  task  at  hand 
since  each  could  perform  the  approximate  transfor  on  on  a  set  of  points 


with  sufficient  speed  to  compete  with  the  performance  of  the  UN1VAC  1108. 

DIGITAL  EQUIPMENT  Corporation's  PDP11  systems  were  not  studied  as  thoroughly 

as  the  other  three  systems,  however,  because  it  was  found  that  they  could 

not  meet  the  timing  specifications  if  rigorous  transformation  was  to  be 

performed  on  more  than  a  small  percentage  (=1%)  of  a  given  set  of  points. 

Based  on  manufacturer's  literature  the  arithmetic  logic  is  much  slower  than 

that  of  the  other  systems  studied,  but  they  can  still  perform  to  specs 

provided  that  few  rigorous  transformations  are  required  (It  should  be  noted 

that  the  current  MARS  system  rigorously  transforms  only  about  0.5%  of  the 
set  of  points,  the  remainder  being  approximately  transformed).  The  sections 

following  provide  a  closer  look  at  the  DATA  GENERAL  Eclipse  S/250  and  the 

PERKIN-ELMER  3220  and  3240  minicomputers.  Conclusions  and  recommendations  are 

given  in  Section  5. 


3. 1  The  DATA  GENERAL  Eclipse  S/250 

The  DATA  GENERAL  Eclipse  S/250  is  a  general  purpose  16  bit  minicomputer 
with  a  large  and  varied  instruction  set.  Its  optional  floating  point  processor 
has  the  capability  of  64  bit  precision,  necessary  for  the  transformations 
used  and  can  preserve  least  significant  digit  accuracy  over  many  compu¬ 
tations.  The  S/250  also  has  a  Writable  Control  Store  option  which  could  prove 
to  be  highly  useful  for  improving  the  speed  of  certain  time  intensive  com¬ 
putations.  With  this  option  installed,  frequently  used  functions  or  subroutines 
(such  as  the  trigonometric  functions  SIN,  COS,  and  TAN)  could  be  coded  to  work 
much  faster  than  that  possible  with  assembly  language.  The  resulting  increase 
in  execution  speed  could  be  significant  (see  Section  4).  In  general,  the  ability 
to  microcode  a  computer  allows  one  to  tailor  the  hardware  to  specific  needs. 

DATA  GENERAL  manufactures  its  own  disk  drives.  The  drives  considered 
have  a  capacity  of  about  50M  bytes  broken  down  into  5  surfaces  each  with  815 
tracks.  Each  track  contains  24  sectors  of  512  bytes.  Though  the  desired  disk 
addressing  would  be  in  powers  of  two,  it  is  not  felt  that  the  effect  of  non¬ 
power  of  two  addressing  present  on  these  drives  will  be  significant.  These 
disk  drives  are  high  speed  with  an  average  access  time  of  43  milliseconds. 

The  DATA  GENERAL  tape  drives  are  industry  compatible  9  track  800/1600 
bpi  75  inch/second,  capable  of  working  with  the  high  speed  Burst  Multiplexor 
Channel  option  of  the  S/250.  The  Burst  Multiplexor  Channel  permits  very 
fast  direct  memory  access  for  both  the  tape  and  disk  drives.  Due  to  its 
relative  low  cost  (=2%  of  total  system  cost)  it  is  felt  that  this  option 
should  be  included  if  the  S/250  is  selected  fo  the  MiniMARS  project. 

One  interesting  aspect  of  the  S/250  is  the  optional  built-in  array 
processor.  The  regridding  and  transformation  task  lends  itself  to  imple¬ 
mentation  on  an  array  processor  since  the  transformation  is  to  be  performed 
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on  a  set  of  points,  not  unlike  operating  with  a  vector  of  elevation  values. 
This  processor,  installed  in  the  main  CPU  cabinet,  will  recognize  special 
instructions  which  allow  the  same  operation  to  performed  on  a  set  of  values. 
The  possible  increase  in  execution  speed  is  similar  to  that  obtained  with 
the  Writable  Control  Store.  The  array  processor  permits  a  floating  point 
precision  of  32  bits,  adequate  for  many  parts  of  the  regridding  and  mosaicking 


3.2  The  PERKIN-ELMER  3220  and  3240 

The  PERKIN-ELMER  3220  minicomputer  and  its  upgrade,  the  3240,  are 
science  oriented  32  bit  computers.  Both  machines  offer  a  floating  point 
processor  option  with  64  bit  precision.  The  major  difference  between  the 
PERKIN-ELMER  floating  point  processor  and  those  of  other  manufacturers 
(including  DATA  GENERAL)  is  the  rounding  method  used  after  performing  a 
floating  point  operation.  Called  R*  (R-STAR)  rounding,  the  method  can 
guarantee  results  at  least  as  accurate  as  those  obtained  by  truncation  or 
one  way  rounding  over  long  series  of  computations.  In  actual  fact,  R*  results 
are  usually  more  accurate  and  could  prove  useful  for  the  transformation 
software. 

Both  the  3220  and  3240  have  multiway  memory  interleaving  which  can 
increase  the  apparent  memory  speed  with  a  resultant  decrease  in  execution 
time.  The  PERKIN-ELMER  3240  has  a  Writable  Control  Store  option  which  can 
cause  significant  improvement  in  execution  speed,  to  the  extent  that  transfor¬ 
mation  times  on  the  3240  with  microcoded  transcendental  functions  are  com¬ 
parable  to  those  obtained  from  a  UNIVAC  1108  (see  Section  4). 

The  disk  drives  used  by  PERKIN-ELMER  .are  supplied  by  Control  Data  Corp. 
They  are  similar  to  .  ..  DATA  GENERAL  drives  in  layout  but  have  a  slightly 
higher  capacity  (67M  bytes)  and  slightly  faster  access  time  (38ms).  The 
higher  speed  may  contribute  about  a  6%  decrease  in  the  I/O  times  when  compared 
to  those  of  DATA  GENERAL.  The  industry  compatible  tape  drives  which  come  with 
the  PERKIN-ELMER  machines  (as  configured)  are  virtually  identical  to  those 
of  DATA  GENERAL. 

PERKIN-ELMER  minicomputers  do  not  have  a  built-in  array  processor  option. 
Instead,  the  array  processors  offered  by  PERKIN-ELMER  are  in  free  standing 
cabinets  and  are  produced  by  Floating  Point  Systems.  These  processors  come  in  a 
variety  of  configurations  including  models  with  either  32  or  64  bit  precision. 


Certain  parts  of  the  transformation  algorithms  require  a  high  degree  of 
precision.  Should  the  array  processor  option  be  desired,  selecting  one  with 
64  bit  precision  would  allow  more  of  the  transformation  algorithm  to  be 
used  in  a  vector  of  values  at  one  time  than  if  a  32  bit  precision  model  is 
chosen.  The  difference  in  cost  between  the  32  bit  and  64  bit  Floating  Point 
Systems  array  processors  would  have  to  be  taken  into  consideration  before  any 
firm  decision  could  be  made  (FPS  recently  announced  their  64  bit  array  pro¬ 
cessor.  The  32  bit  processor  is  in  the  $100K  range.  The  64  bit  processor  is 
nearly  double  that.). 
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Table  1. 

COMPARATIVE  PRICES 


ITEM 

DATA  GENERAL 
S/250 

PERKIN-ELMER 

3220 

PERKIN-ELMER 

3240 

DLC 
11/2  / 

CPU  +  256K 

Byte  Memory 

$34,425 

$28,440 

$87,690 

$23,545 

Floating 

Point  roc. 

$5,266 

$4,424 

$7,505 

$2635 

High  Speed  (1) 

$18,530 

$12,008 

$12,008 

$21,845 

Disk  Drive  (2) 

$18,530 

$12,008 

$12,008 

$16,405 

Tape  Drive  (1) 

$13,175 

$15,089 

$15,089 

$17,170 

High  Speed  (2) 

$9,605 

$7,603 

$7,603 

$10,880 

Printer 

$2,252 

$3,476 

$3,476 

$1445 

(cabinet) 

Other:  Cache,  $2,600 

„  c  ,  (BMC) 

Power  Supply,  $2>678 

Cabinets,  Clock,  (cache) 

I/O  Interface, etc.  $^>570 

(cabinet) 

$869 

(channels) 

$2,765 

(cache) 

$2,765 

(power) 

$1,343 

(clock) 

$711 

(channels) 

$2,765 

(cache) 

$5,925 

(cabinet 

$1105 

(expansion 

$2720 

(interface 

$731 

(cable) 

$51 

(clock) 

$697 

Software 

$4,400 

$12,600 

$12,600 

$8800 

TOTAL 

$115,031 

$104,101 

$166,669 

$107,978 

*NOTE*  The  PERKIN-ELMER  software  costs  reflect  their  pricing  for  their  very 
powerful  optimizing  FORTRAN  compiler.  Duplication  costs  (for  multiple  systems) 
are  one  tenth  of  initial  costs  (for  example,  $12,600  for  the  first  minicomputer 
system,  $1260  for  each  next  system.) 
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The  costs  of  the  systems  previously  described  are  roughly  equivalent. 

The  16  bit  DATA  GENERAL  Eclipse  S/250,  DEC's  PDP11/44,  and  the  PERKIN-ELMER 
3220  and  3240  are  within  the  $100K  to  $200K  range  mentioned  in  the  initial 
proposal.  The  PERKIN-ELMER  3240,  while  considerably  more  expensive  than  the 
other  systems,  has  high  speed  capabilities  approaching  those  obtainable  with 
a  third  generation  machine  like  the  UNIVAC  1100/81.  All  four  systems  can 
perform  the  proposed  regridding  and  mosaicking  task  within  the  time  constraint 
if  the  approximate  transformation  is  used  for  the  major  portion  of  the  points 
as  is  currently  done  with  the  MARS  system.  The  decision  on  which  machine  to 
choose  should  be  based  not  only  on  the  capability  to  perform  the  regridding 
and  mosaicking  task,  but  also  on  the  capability  to  perform  future  tasks,  the 
ease  with  which  the  system  can  be  integrated  into  the  existing  DMA  framework, 
and  the  reputation  of  the  manufacturer (s)  for  delivery  and  support. 

The  DATA  GENERAL  S/250  is  the  same  basic  machine  used  in  the  UNAMACE 
equipment.  As  such,  DMA  has  personnel  who  are  familiar  with  the  inevitable 
idiosyncrasies  of  a  particular  hardware  system.  However,  the  S/250  is  only 
a  16  bit  machine.  Future  work  which  DMA  m  ay  desire  to  perform  on  the  chosen 
system  could  possibly  not  be  performed  adequately  on  a  16  bit  machine.  On 
the  other  hand,  since  map  oriented  operations  can  quite  often  be  realized 
as  operations  on  matrices  of  elevation  (or  other)  values,  the  built-in 
array  processor  of  the  S/250  could  prove  to  be  very  useful  for  quick  operation* 
on  a  vector  of  values.  The  S/250  is  a  fast  and  powerful  16  bit  minicomputer 
backed  by  a  large  experienced  corporation. 

The  PERKIN-ELMER  3220  (or  3240)  was  designed  as  a  scientific  use  32  bit 
machine.  Its  instruction  set  and  hardware  configuration  are  oriented  towards 


manipulation  of  numerical  data  in  an  efficient  manner.  -In  addition  the  3220 
has  a  32  bit  architecture.  Future  DMA  needs  could  be  implemented  on  this  machine 
with  minimum  modification.  PERKIN-ELMER' s  optimizing  FORTRAN  compiler  produces 
object  code  which  approaches  that  of  an  experienced  programmer.  The  "deluxe" 

3240  minicomputer  is  the  most  powerful  of  the  machines  reviewed  for  this  report. 
Its  extremely  high  speed  operations  combined  with  a  large  memory  address  space 
would  permit  easy  future  software  expansion  according  to  DMA  needs.  PERKIN- 
ELMER  is  a  large  well  known  firm  involved  in  many  different  fields. 

DIGITAL  EQUIPMENT'S  PDP11/44  (or  similar  models)  are  16  bit  minicomputers 
with  general  purpose  instruction  sets.  The  speed  of  this  machine  is  adequate 
when  an  approximate  transformation  is  used.  Being  one  of  the  first 
minicomputers  introduced, the  PDP  line  has  a  proven  architecture  for  a  variety 
of  tasks.  In  addition,  the  PDP  is  familiar  to  more  programmers  than  the  other 
machines,  largely  due  to  its  presence  at  many  universities.  DEC  is  the  largest 
minicomputer  manufacturer  and  is  known  for  its  support  services. 

As  this  report  was  in  the  final  stages  of  being  written  DEC  announced 
a  less  expensive  version  of  their  VAX  line  of  computers.  The  VAX  is  an 
extremely  powerful  32  bit  machine  with  a  large  memory  address  space.  Before 
any  decision  is  made  on  the  type  of  machine  to  be  used  for  the  MiniMARS 


system  we  feel  that  this  machine  should  be  examined. 
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4 .  BENCHMARKS 


The  MiniMARS  algorithm  relies  on  the  ability  of  a  minicomputer  system 
to  perform  certain  operations  such  that  when  comparisons  are  made  the  time 
required  on  a  UNIVAC  1100/81  and  the  minicomputer  are  not  radically  different. 
With  this  requirement  in  mind,  benchmark  tests  were  performed  on  a  UNIVAC 
and  on  four  different  minicomputers.  Table  2  shows  execution  times  for 
three  different  transformations  performed  on  10^  points  and  the  extrapolated 
time  estimates  for  10^  points.  The  tests  were  performed  on  a  UNIVAC  1108 
and  on  a  UNIVAC  1180.  The  minimal  driver  and  transformation  routines  were 
then  transported  to  a  PDP11/45,  P/E  3220  and  3240,  and  DG  S/250  machines. 

Table  3  gives  the  times  for  10^  points  and  extrapolated  values  for  10^ 
points.  As  can  be  seen,  the  times  for  the  obsolete  PDP11/45  are  inadequate 
for  rigorous  transformation  of  all  points  of  a  model. 

The  PERKIN-ELMER  machines  show  two  times.  These  result  from  using  a 
development  compiler  and  an  optimizing  compiler.  The  times  for  the  3220 
are  within  a  factor  of  three  of  those  of  the  UNIVAC  1100/81.  The  times  for 
the  3240  are  within  a  factor  of  two  of  those  obtained  from  the  UNIVAC  sys¬ 
tems,  demonstrating  the  speed  with  which  state-of-the-art  minicomputers  can 
perform.  The  lower  estimate  of  the  3240  machine  was  obtained  by  using  a 
compiler  which  places  certain  operations  in  writeable  control  store  (in 
this  case,  sines  and  cosines).  The  use  of  microcoded  instructions  can 
increase  performance  by  as  much  as  30%. 

In  order  to  test  I/O  speeds  for  disk  accessing,  a  test  was  run  on  a 
PDP11/45  running  Bell  Laboratory's  UNIX  operating  system.  The  disk  used 
was  an  RK05  equivalent  with  70  ms  average  ac c  is  time.  This  increase  in 
access  time  over  the  proposed  access  time  should  be  taken  into  considera¬ 
tion.  Table  4  gives  the  results  of  4  tests  of  pseudo-random  reads  of 


4-1 


Table  2. 

Transformation  Benchmarks  For  MiniMARS 


UNIVAC 

1108 

UNIVAC  1 

100/81 

Transformation 

CPU  Time  for 

4 

10  Points 

Extrapolated 
10^  Points 

CPU  Time  for 

4 

10  Points 

Extrapolated  ; 

6  1 

10  Points  | 

Geographic  to 

UTM 

7.1  seconds 

12  minutes 

8.5  seconds 

1 

1 

14  minutes  j 

UNAMACE  to 

Geographic 

13.8  seconds 

23  minutes 

16.5  seconds 

t 

28  minutes  j 

UNAMACE  to 

UTM 

20.1  seconds 

35  minutes 

24.6  seconds 

j 

41  minutes 

j 

*  —  These  time  were  obtained  while  running  the  tests  in  "DEMAND"  mode  on 
the  machine  at  DMAHTC.  The  algorithm  used  for  computing  execution  time  has  i 

been  improved  on  the  1100/81.  This  probably  accounts  for  the  longer  execution 
times  than  on  the  UNI VAC  1108. 

Table  3. 

Geographic  to  UTM  Times  on  Four  Minicomputers 


DATA 

GENERAL 

S/250 

PERKIN  * 

ELMER 

3220 

PERK] 

ELMI 

324C 

* 

IN 

HR 

> 

DEC  | 

PDP  i 

11/45  ! 

CPU  Time  for 

4 

10  Points 
(seconds) 

20 

30 

23 

17.6 

t 

12.8 

132 

Extrapolated 
for  10^  Points 

(minutes) 

_ 

33.3 

50 

36.3 

1 

t 

— I 

28 

! 

21.3 

220  | 

1 

I 

Two  times  are  shown  for  the  Perkin-Elmer  machines.  The  longer  time  in  each 
case  was  obtained  using  a  "development"  compiler,  the  shorter  time  was  obtained 
using  an  optimizing  compiler. 

The  analysis  in  the  original  proposal  used  a  range  of  3  to  45  minutes  for 
10^  points  to  be  transformed.  If  rigorous  ins f ormat ion  is  going  to  be  used  for 

each  point  only  the  I’DPll/43  could  not  perform  adequately. 
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10^  512  byte  blocks.  Table  5  contains  the  results  of  4  tests  of  10^*  sequen¬ 
tial  reads  of  512  contiguous  blocks.  Projecting  the  timing  values  obtained 
into  those  possible  with  the  specified  high  speed  disk  demonstrate  that  the 
I/O  operations  can  be  performed  with  sufficient  speed,  especially  if  there 
is  concurrent  input  and  output. 

The  benchmark  tests  have  shown  that  current  minicomputers  do  have  the 
power  and  speed  to  compete  with  large  main  frames  such  as  the  UNIVAC  1100/81 
The  high  reliability  and  low  cost  of  minicomputers  strongly  suggest  that 
they  should  be  considered  not  only  for  the  implementation  of  the  algorithm 
in  this  report,  but  also  for  other  algorithms  performed  by  the  MARS  system. 

the  implementation  of  this 
specialized  hardware. 


The  following  section  gives  recommendations  on 
and  other  algorithms  on  both  minicomputers  and 
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5 .  CONCLUSIONS  AND  RECOMMENDATIONS 
5.1  General  Conclusions 

(Gl)  The  MARS  regridding  step  can  be  performed  in  a  very 

simple  manner  as  originally  proposed,  at  a  rate  of  10^  points 
in  40  to  80  minutes,  and  on  hardware  costing  roughly  $125,000. 

This  conclusion  is  based  on  the  decomposition  of  the  processing  into 
simple  steps  each  one  of  which  is  tractable  to  analysis  and  can  be  cali¬ 
brated  by  benchmarks.  Actual  benchmarks  for  I/O  and  computation,  reported 
in  Section  4,  support  the  previous  conclusion. 

(G2)  The  current  MARS  software  system  could  be  transported  to  a 

minicomputer  system  with  no  loss  in  execution  speed  and  a  great 
saving  in  execution  cost. 

The  minicomputers  studied  can  have  large  memories  and  fast  arithmetic 
and  could  easily  host  the  current  MARS.  Throughput  should  increase  if  the 
system  were  dedicated  to  MARS  through  decreased  overhead  and  I/O  contention 

Moreover,  memory  limitations  are  known  to  slow  down  the  current  UNIVAC  MARS 
Large  memory  increases  possible  in  the  mini  systems  could  speed  up  the  pro¬ 
cessing  by  reducing  fragmentation.  Costs  of  converting  the  current  UNIVAC 
system  would  be  high  because  of  the  complexity  of  that  system  -  probably 
higher  than  the  cost  of  creating  the  proposed  Mini-MARS  software  system. 
Conversion  of  the  current  MARS  was  not  actually  a  task  of  the  current  pro¬ 
ject,  however,  given  the  current  information  it  is  easy  to  evaluate  such 
an  option. 
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5.2  Conclusions  Regarding  Coordinate  Transformations 

(Tl)  Rigorous  transformation  of  every  grid  point  can  double  execution 
time  relative  to  regridding  with  a  large  amount  of  interpolation. 

This  conclusion  is  based  on  data  from  the  MARS  system  contained  in  the 
work  proposal  and  confirmed  by  subsequent  analysis  and  benchmarks.  Actually, 
rigorous  transformation  creates  a  50/50  mix  of  computation  and  I/O  time. 
Minicomputers  with  slower  arithmetic  but  comparably  fast  I/O  would  cause  a 
strong  shift  toward  being  compute  bound.  See  Appendix  A  for  a  study  of  the 
MARS  transformation  routines. 

(T2)  Rigorous  transformation  should  be  avoided  on  all  systems 
considered. 

From  Tl  it  is  seen  that  even  on  very  fast  conventional  computers  (i.e. 
ITNIVAC  1108,  PERKIN-ELMER  3240)  half  of  the  process  time  will  be  taken  up  in 
transformation.  Of  the  computers  studied  only  the  UNIVAC  could  meet  the  80 
minute  process  time  if  all  points  were  regridded  rigorously. 

The  PDP  ll's  would  require  several  hours  for  rigorous  transformation  but 
could  finish  in  less  than  80  minutes  using  interpolation. 

(T3)  Provided  that  only  1  of  200  grid  points  were  transformed  rigorously 
all  of  the  computers  studied  (PDP  11,  S/250,  UNIVAC  1108,  3220  and 
3240)  could  regrid  10^  points  in  less  than  80  minutes. 

This  is  a  direct  consequence  of  the  analysis  of  the  algorithms  and  the 
timings  from  Tables  2  and  3.  Past  experience  with  MARS  indicates  that  the 
fraction  1  of  200  is  very  reason  le;  we  will  seldom  have  to  rigorously  trans¬ 
form  more  points  than  this. 

(T4)  At  this  point  in  time  special  off-the-shelf  hardware  to  speed 
up  the  transformations  does  not  appear  to  be  cost  effective. 
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The  DG  array  processor,  in  the  $14,000  range,  is  the  only  inexpensive 
array  processor  available  for  any  of  the  studied  machines  but  it  does  not 
handle  double  precision  and  is  not  a  true  parallel  processor.  A  recently 
announced  Floating  Point  64-bit  array  processor  is  in  the  $150,000  - 
$200,000  range.  Since  such  a  device  could  only  half  the  computation  time, 
for  example  on  the  PE3240,  it  would  be  wiser  to  double  capacity  by  having 
two  conventional  systems.  Future  reductions  in  cost  and  interfacing  pro¬ 
blems  would  cause  us  to  reconsider. 

(T5)  The  costs  of  the  16  bit  and  32  bit  machines  are  close  enough 
together  that  the  selection  of  the  32  bit  machine  appears  to 
be  best.  (This  machine  can  be  either  the  3220  or  3240,  or 
one  not  reviewed  by  L.N.K.). 

The  capability  for  future  expansion  and  the  easy  addressing  methods  of 
the  32  bit  CPU  are  very  strong  assets. 


5.3  Final  Conclusions 


The  Mini-MARS  project  can  be  implemented  on  a  mini-computer .  It  is 
recommended  that  a  32  bit  machine  be  selected.  The  MARS  system  even  in  its 
present  form  could  be  implemented  on  the  same  machine  (as  it  now  exists)  with 
some  increase  in  thru-put.  Hardware  costs  will  be  in  the  proposed  range  of 
$100K  to  $200K.  Setup  and  software  costs  should  be  within  S100K  to  S200K. 
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STUDY  OF  DMA  TRANSFORMATION  PACKAGE 

The  DMA  coordinate  transformation  package  contains  several  superb  algorithms 
for  performing  coordinate  conversion.  Two  subroutines,  however,  require  modi¬ 
fication  in  order  to  avoid  possible  incorrect  computation.  The  suggested  changes 
are  given  in  Section  A. 

If  a  large  set  of  points  (10,000  points,  for  instance)  are  to  be  processed 
using  the  subroutines,  some  reduction  in  computation  time  can  be  achieved  by 
avoiding  recomputing  constants  such  as  n / 2 ,  the  eccentricity  of  the  spheroid, 
coefficients  for  equations ,  etc.,  in  each  call  of  a  subroutine.  The  number  of 
exponentiations,  multiplications  and  divisions  should  be  minimized.  Specific 
suggested  modifications  in  J:his  direction  are  given  in  Section  B.  It  should  be 
pointed  out  that  the  saving  in  computation  time  using  such  modification  is  not 
likely  to  exceed  10%.  On  the  otherhand,  if  such  economization  is  combined  with 
improvements  in  the  computation  of  tht  trigonometric  functions,  the  total  saving 
in  computing  time  can  exceed  25%.  In  Section  C,  methods  for  reducing  the  compu¬ 
tation  time  of  the  trigonometric  functions  are  discussed.  The  actual  FORTRAN 
package  is  given  at  the  end  of  this  memo. 

A.  Corrections  and  Safeguard  Steps 

1.  FLHUPS,  line  9:  Change  "CON  =  1.0"  to  "CON  =  2.0D0". 

The  derivation  of  the  value  for  CON  can  be  found  in  Adter  (  [1]  ),  page  97. 

2.  FLHUPS :  Change  all  occurences  of  the  FORTRAN  library  functions  ABS,  SIN, 
COS  and  SQRT,  to  DABS,  DSIN,  DCOS  AND  DSQRT,  respectively.  This  is  a 
precautionary  step. 

3.  FLHXYZ,  line  3:  Change  "C  THIS  PROGRAM  TRANSFORMS  GEOGRAPHIC  COORDINATES 
TO  LOCAL"  to  "C  THIS  PROGRAM  TRANSFORMS  GEOGRAPHIC  COORDINATES  TO 


GEOCENTRIC". 


4.  UPSFLH,  line  10:  Change  "CON  =  1.0"  to  "CON  =  2.0D0". 

5.  UPSFLH:  Clange  all  occurences  of  the  FORTRAN  library  functions  ATAN,  ABS, 

SIN,  COS  and  SQRT  to  DATAN,  DABS,  DSIN,  DCOS  and  DSQRT,  respectively. 

5  XYZFLH,  line  28:  The  variable  "A"  should  be  assigned  the  correct  value 

before  this  step  is  executed.  Such  an  assignment  step  seems  to  be  missing. 
Also  the  logic  in  the  expression  itself 

FLH(3)  =  A  -  ABS (Z) 

appears  to  contradict  the  definition  of  geographic  coordinates  (unless  this 
is  a  special  convention  used  at  DMA).  The  formula  one  would  generally  use  is 

FLH(3)  =  DABS (Z)  -  AXES(2). 

7.  UTMFLH:  The  Newton  i&ration  loop  (lines  33-37)  should  contain  an  accuracy 
check  in  order  to  (i)  alert  the  user  to  possible  inaccurate  results,  and 
(ii)  eliminate  unnecessary  loops  when  an  accurate  result  has  been  obtained. 
See  lines  47-50  of  UPSFLH  for  an  example  of  such  a  checking  device. 


B .  Reducing  Computation  Time  (Economization  of  Arithmetric  Operations) 

1.  In  an  ini t ial iza t ion  subroutine,  initialize  PI  and  PIHALF: 

PI  =  3. 14 159265 3589 79 3D0 
PIHALF  =  PI/2.0D0 

Save  the  results  in  a  COMMON  area. 

A 


Change  every  occurence  of  "PI/2.0"  to  "PIHALF": 

line  14  of  FLHUPS , 

line  51  of  UPSFLH, 

lines  21,  25,  27,  30  of  XYZFLH. 

In  an  initialization  subroutine,  compute  the  eccentricity  E  and  its 
square  E2: 

C  A  is  semi-major  axis  and  B  is  semi-minor  axis 

A  =  . 

B  -  . 

E2  =  (  A  +  B  )*(  A  -  B  )/A/A 
E  =  DSQRT(E2) 

Save  the  values  of  E  and  E2  in  a  COMMON  area. 

Remove  lines  12,  13  of  FLHUPS, 

line  30  of  FLHUTM, 
lines  12,  13  of  FLHXYZ, 
lines  25,  26  of  UPSFLH, 
line  20  of  UTMFLH, 
line  41  of  XYZFLH. 

In  an  initialization  subroutine,  compute  the  coefficients  A,  B,  C  and  D 
(lines  32-35  of  FLHUTM) .  Save  the  results  in  a  COMMON  area.  E  and  E2 
should  be  computed  before  computing  A,  B,  C  and  D. 

Change  lines  19-21  of  FLHUPS  to 

TEMPL  -  (1.0D0  -  E)/ (1.0D0  +  E) 

R  -  CON*(A/B)*TANZ*A*(TEMP*TEMPL)**EX. 
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Change  line  43  of  UPSFLH  to 

TEMP 2  =  TEMPI/TEMP. 

In  each  case,  an  exponentiation  operation  is  eliminated.  Exponentiation 
is  costly:  approximately  eight  multiplications  are  necessary  to  perform 
an  exponentiation  operation. 

5.  In  UPSFLH,  replace  lines  34  and  35  with 

S2  =  2.0D0  *  SI  *  Cl 

C2  =  2.0D0  *  Cl  *  Cl  -  1.0D0. 

This  modification  eliminates  one  call  of  the  function  DSIN(P)  and  one 
call  of  the  function  DCOS(P). 

6.  In  XYZFLH,  add  the  line  "RXY  =  X*X  +  Y*Y"  after  the  line  42. 

Change  line  45  to  "R  =  RXY  +  ZP*ZP",  and  change  line  54  to  "R  =  RXY'" . 

C.  Computation  of  Trigonometric  Functions 

In  many  subroutines,  the  time  required  to  compute  the  trigonometric 
functions  is  significant.  As  a  measure  of  this  segment  of  computing  time, 
we  estimate  the  ratio  R: 

R  =  (Total  number  of  multiplications  and  divisions  necessary  to  compute 
the  trigonometric  functions  in  subroutine)/ (Total  number  of  multi¬ 
plications  and  divisions  performed  in  each  call  of  the  subroutine) 

If  a  DO-loop  is  to  be  executed  in  the  subroutine,  the  ratio  R  takes  into 
account  only  the  operation  counts  involved  in  one  loop.  The  ratio  will 
decrease  somewhat  if  more  than  one  loop  is  considered.  The  economization 
modifications  of  Section  B  are  assumed. 
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Fortran  trigonometric  functions  are  usually  implemented  using  Chebvshev 
polynomials  or  Lagrange  interpolation  polynomials.  An  example  of  an  approxi¬ 


mating  polynomial  for  cos  ^ttx  is 

.9999999724  -  .3084242527  x2  +  .0158499130  x4  -  .0003188790  x6  , 

Fike  (  [2]  ),  page  133.  In  this  case,  four  multiplications  are  performed 
during  the  evaluation  of  the  polynomial  using  the  factored  form 
a  +  (  (b  +  (c  +  dx^)x^  )x^  with  a,  b,  c,  d  as  above. 


SUBROUTINE  R  (approximate) 

FLHUPS  -  50% 

FLHUTM  -  33% 

FLHXYZ  -  50% 

LOGO  0% 

(LOGO  calls  FLHXYZ  and  0RTM6A) 

MOVE  0% 

(MOVE  is  called  by  TRANCD  only) 

TRANCD  -  0% 


(TRANCD  is  an  administrative  subroutine.  It  calls  other 
subroutines. ) 


UPSFLH  -  35% 

UTMFLH  -  40% 

XYZFLH  -  20% 


The  above  ratios  justify  considering  "hard-wiring"  the  trigonometric  functions, 
implementing  faster  piecewise  polynomial  approximations. 

Note  the  following  computation  times  cited  for  DSIN  along  with  the  source. 

25  x  10  ^  second  DG  S/250  Eclipse  with  microcoded  functions.  Timing 


gotten  from  technical  literature. 
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40  x  10  second  Univac  1108.  Timing  gotten  from  actual  benchmark 
program  run. 

97  x  10  k  second  Perkin  Elmer  3220.  Timing  gotten  from  technical 
literature. 

4 

Benchmarks  actually  run  with  10  points  transformed  from  geographic  to 
UTM  coordinates  are  as  follows. 

7.]  sec  Univac  1108 

23.0  sec  Perkin  Elmer  3220 

12.8  sec  Perkin  Elmer  3240 

D.  Documentation 

Specific  source-references  should  be  given  for  the  approximation  formulas 
in  FLHUTM  and  UTMFLH,  since  the  formulas  in  the  two  subroutines  are  modifi¬ 
cations  of  the  original  defining  equations  of  the  UTM  coordinates. 
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DMA  TRANSFORMATION  SOFTWARE 

Following  is  a  listing  of  all  transformation  software  currently  in  use  by 
DMA.  Line  numbers  in  these  listings  are  those  referenced  in  the  prior  text. 

Four  driver  programs  are  listed  before  the  subroutines  in  the  package.  These 
drivers,  perhaps  with  some  changes  in  constants,  were  used  to  run  benchmarks  on 
the  transformation  software.  Benchmarks  were  run  on  two  different  Univac 
installations  (1110  and  118t)  and  two  different  Perkin  Elmer  installations 
(3220  and  3240). 
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000001 

000 

c 

DRIVER  TO  CAL*  TRAMCB  FOR  6E06RAPHIC  TO  UTN  CONVERSION 

000002 

000 

c 

THIS  DRIVER  SHORT  CUTS  BY  NOT  CALLING  ROUTINE  TRANCD***** 

000003 

000 

c 

ALSO  WRITTEN  'UNIVAC  FREE'  TO  GO  TO  PERK1N-ELHER  FOR  BENCHMARK 

000004 

000 

c 

000005 

000 

IMPLICIT  INTEGER ( A-Z > 

000006 

000 

DOUBLE  PRECISION  PTIN(3>,PT0UT(3),DAXES(2),P0(3), HO, FLAG, ZONE 

000007 

000 

DOUBLE  PRECISION  F(4> 

000008 

000 

c 

000009 

000 

MOUNT  =  0 

000010 

000 

DPI  *  3. 1 4 1592653589793 24 DO 

000011 

000 

PTIN(I)  «  10. *DP 1/180. 

000012 

000 

PT I N ( 2 )  «  1 0.*DPI/1 80. 

000013 

000 

PTINI3)  -  O.ODO 

000014 

000 

BAXES ( 1 )  >  637813SD0 

000015 

000 

DAXES ( 2 )  =  6356750. 52D0 

000016 

000 

c 

CODE  ADDED  TO  SHUNT  C.UL  TO  TRANCD 

000017 

000 

F<3)=DAXES<1> 

000018 

000 

F(4)*DAXES(2) 

000019 

000 

pom  «  io.odo 

000020 

000 

P0( 2 )  =  20.*DPI/180. 

000021 

000 

P0<3)  =  100.*0PI/180. 

000022 

000 

HO  *  O.ODO 

000023 

000 

FLAG  *  O.ODO 

000024 

000 

ZONE  «  O.ODO 

000025 

000 

KFLAG  *  7 

000026 

000 

DO  50  I  =  1,200 

000027 

000 

DO  50  J*1 ,50 

000028 

000 

PTIN< 1 )  *  PT I N < 1 )  +  1 . *DP 1/180. 

000029 

000 

PTIN(2)  *  PTIN(2)*t.*DPI/180. 

000030 

000 

PT IN ( 3 )  =  PTIN(3)*1 . 

000031 

000 

KOUNT  «  MOUNT  +  1 

000032 

000 

c 

CALL  TRANCD (PTIN, PTOUT, DAXES, PO, HO, FLAG, ZONE, MFLAG) 

000033 

000 

c 

HERE  IS  CODE  DONE  IN  TRANCD  ENABLING  THAT  ROUTINE  TO  BE  SHUNTED 

000034 

000 

F(1)*PTIN(1) 

000035 

000 

F(2)*PTIN<2) 

000036 

000 

CALL  FLHUTM1F, PTOUT, ZONE) 

000037 

000 

PTOUT  <3)*PTIN<3) 

000038 

000 

IF  ( I . Efl . 1 . AND. J.EO. 1 )  UR I TE (6,49)  PTIN, PTOUT 

000039 

000 

IF( I .EQ. 200. AND. J.EQ .50 >  UR I TE < 6 , 49 )  PTIN, PTOUT 

000040 

000 

49  FORMAT < '  PTIH  '.3D20.10,/,'  PTOUT  ',3D20.10,//> 

000041 

000 

5C 

1  CONTINUE 

000042 

000 

URITE46, 100)M0UNT 

000043 

000 

100  FORMA T< 'ONORMAL  EXIT.  COUNT  *  ',16) 

000044 

000 

STOP 

000045 

000 

END 
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1  C  DRIVER  TO  CALL  TRANCD  FOR  UNAHACE  TO  GEOGRAPHIC  CONVERSION 


2 

C 

3 

IMPLICIT  INTEGER(A-Z) 

4 

DOUBLE  PRECISION  PTIN<3) ,PT0UT(3),DAXES<2> , POO), HO, FLAG, ZONE, DPI 

5 

C 

A 

KOUNT  =  0 

; 

PTIN(I)  =  -5000. ODO 

8 

PTIN12)  '  -5000. ODO 

V 

DAXES11  )  =  6378135D0 

10 

DAXES  <  2 )  =  6356750. 52D0 

ii 

DPI  =  3. 1 4 1 5926535897 9 324D0 

12 

POd)  =  0.0  ^LATITUDE 

13 

PO  <  2 )  =  20 .  *DP 1/1 80  .  8L0NGI TUDE 

14 

P  0  C  3 )  =  1 00. *DPI/1 80.  8 AZINUTH 

15 

HO  =  0. ODO 

16 

FLAG  =  O.ODO 

17 

ZONE  =  O.ODO 

18 

KFLAG  =  16  0  FOR  UNA  TO  GEO 

19 

DO  50  I  =  1,10 

20 

P  T IN ( 3 )  =  O.ODO 

21 

POd  )  =  POd )+5. «DPI / 1 80 . 

22 

DO  50  J  =  1,10000 

23 

PTINd  )  =  PTINd  )  +  1. 

24 

PTINI2)  =  PTIN1 2  )  +  1 . 

25 

PTIN ( 3 )  =  PTIN(3)+1 . 

26 

KOUNT  =  KOUNT  ♦  1 

27 

50 

CALL  TRANCD (PTIN, PTOUT, DAXES, PO, HO, FLAG, ZONE, KFLAG) 

28 

URI TE ( 6 , 1 00 ) KOUNT 

29 

100 

FORMAT ( 'ONORMAL  EXIT,  COUNT  =  ',16) 

30 

STOP 

31 

END 
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2 

C 

3 

IMPLICIT  INTEGER < A-Z > 

4 

DOUBLE  PRECISION  PTIN(3> ,PT0UT (3 ) ,DAXES(2) ,P0(3) , HO, FLAG, 

5 

C 

4 

KOUNT  =  0 

7 

PTINd)  *  -5000. ODO 

8 

PTIN12)  =  -5000. ODO 

9 

DAXESd)  *  437813SD0 

to 

DAXES12)  =  4354730. 32D0 

11 

DPI  *  3. 141 592653589 79324D0 

12 

POO)  *  0.0  PLATITUDE 

13 

P0(2)  =  20.*DPI/1 80.  8L0NGITUDE 

14 

P0<3)  =  1 00.*DPI/1 80 .  8AZIHUTH 

15 

HO  =  O.ODO 

16 

FLAG  =  O.ODO 

17 

ZONE  =  O.ODO 

18 

KFLAG  =17  P  FOR  UNA  TO  UTN 

1? 

DO  50  I  =  1,10 

20 

PTIN13)  =  O.ODO 

21 

P0(1)  =  POd  )+5. ♦DPI/180. 

22 

DO  50  J  =  1,10000 

23 

PTINd)  =  PT I N  <  1  )  +  1 . 

24 

PTIN12)  *  PTIN(2)  +  1 . 

23 

PTIn  <  3 )  =  PTIN(3)  +  1 . 

26 

KOUNT  =  KOUNT  ♦  1  17 

27 

50 

CALL  TRANCD (PTIN, PTOUT,DAXES,PO, HO, FLAG, 20NE, KFLAG) 

28 

URITE (4, 1 001K0UNT 

29 

100 

FORNAT ( 'ONQRHAL  EXIT,  COUNT  *  M6> 

30 

STOP 

31 

END 

1 

C  DRIVER  TO  CALL  TRANCD  FOR  GEOGRAPHIC  TO  UTH  CONVERSION 

2 

C 

3 

IMPLICIT  INTEGER! A-Z ) 

4 

DOUBLE  PRECISION  PT IN < 3 ) ,PTOUT ( 3 ) , DAXES < 2 ) ,P0(3) , HO, FLAG 

5 

C 

6 

KOUNT  =  0 

7 

/ 

'Pi  IN(  1 )  =  10. ♦DPI/180. 

8 

' 

PTIN(2)  =  10. ♦DPI/180. 

9 

1 

PTIN13)  =  O.ODO 

10 

DAXES ( 1 )  =  6378135D0 

11 

v 

QAXES ( 2 )  =  6354730. 52D0 

12 

^DPI  *  3. 1 4139245358979324D0 

13 

POd)  *  10. ODO  PLATITUDE 

14 

PO ( 2 )  <  20. ♦DPI/180.  PLONGITUDE 

13 

POO)  *  100. ♦DPI/180.  8AZIHUTH 

16 

HO  =  O.ODO 

17 

FLAG  =  O.ODO 

18 

ZONE  =  O.ODO 

19 

KFLAG  =  7  P  FOR  6E0  TO  I"* 

20 

DO  50  I  =  1,2000 

21 

DO  50  J  =  1,50 

22 

PTINd)  *  PT1N<1)+1. ♦DPI/180. 

23 

PT  I N  <  2 )  =  PTIN<2>+1. ♦DPI/180. 

24 

PT I N ( 3 )  «  PT I N <  3 ) ♦ 1 . 

25 

KOUNT  *  KOUNT  ♦  1 

26 

50 

CALL  TRANCD  (PTIN, PTOUT, DAXES, PO, HO, FLAG, ZONE,  KFLAG) 

27 

URITE  <6, 1 00 )KOUNT 

28 

100 

FORNAT ( 'ONQRHAL  EXIT,  COUNT  =  ',16) 

29 

STOP 

30 

END 

A-  1  1 


r 


I 

SUBROUTINE  FLHUPS  <PTIN,PTOUT) 

FLH00100 

2 

C 

FLH00200 

3 

C 

THIS  ROUTINE  GEOGRAPHIC  COORDINATES  TO  UNIVERSAL  POLAR  STEREOGRAPHIC 

FLH00300 

4 

C 

FLH00400 

3 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

FLH00500 

6 

DIMENSION  PTIN<5) , PTOU T (2) 

FLH00600 

7 

C 

FLH00700 

e 

PI  ■  3.141592653589793 

FLH00800 

9 

CON  =  1.0 

FLH00900 

10 

A  =  PTIN<4) 

FLH01000 

11 

B  -  PTIN(5) 

FLH01100 

12 

E  *  1 .  -  (B/A)**2 

FLH01200 

13 

E  *  SORT  (E) 

FLH01300 

H 

P  *  PI/2.0  -  ABS  (PTINd)) 

FLH01400 

13 

TANZ  *  SIN  (P/2.0)  /  COS  (P/2.0) 

FLH01500 

16 

COSP  =  COS  (P) 

FLN01600 

17 

TEMP  =  <1.  +  £  *  COSP)  /  <1.  -  E  ♦  COSP) 

FLH01700 

18 

EX  =  E/2. 

FLH01800 

19 

TANZ  *  TANZ  *  TEMP**EX 

FLH01900 

20 

TEMP  =  (1.  -  E)  /  (1.  ♦  E) 

FLH02000 

21 

R  *  CON  *  (A/B)  *  TANZ  *  A  «  TEMP**EX 

FLH02100 

22 

PTOUT < 1 )  =  R  •  COS  <PTIN(2>)  ♦  2000000.0 

FLH02200 

23 

Y  =  R  *  SIN  <PTIN<2) ) 

FLH02300 

24 

IF  (PTIN(D)  10,20,20 

FLH02400 

25 

10 

PTOUT (2)  =  2000000.0  ♦  Y 

FLH02500 

26 

GO  TO  30 

FLH02600 

27 

20 

PTOUT (2)  =  2000000.0  -  Y 

FLH02700 

28 

30 

RETURN 

FLH02800 

29 

END 

FLH02900 
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a...  j 


FLHUTM  1  of  2 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 
13 
16 

17 

18 
1? 


SUBROUTINE  FLHUTHIF, 0.Z0NE ) 

C  6E06RAPHIC  TO  UTH 

IMPLICIT  DOUBLE  PRECISION  <A-H,L-Z) 

DIMENSION  F(4),  0(2) 

C  F  AND  0  ARE  THE  SAME  CALLINS  VARIABLES  THAT  APPEAR  IN  GENCOR 
C  F< 1 )  (OR  LI)  IS  THE  LATITUDE  IN  RADIANS 

C  F < 2 )  (OR  L2)  IS  THE  LONGITUDE  IN  RADIANS 

C  F ( 3)  (OR  AD  AND  F(4)  (OR  B1 >  ARE  THE  AXES  OF  THE  SPHEROID. 

LI *F ( 1 ) 

L2*F (2) 

A1*F(3) 

I1«F(4) 

Z  «  ZONE 

C  Y  IS  A  CONVERSION  FACTOR  DEG  TO  RAD 
Y *.01 743329251 99 4 32958D0 
IF  (ZONE .6T.60.0D0)  ZONE  *  O.ODO 
L0N6  *  L2/Y 

IF  (LONG. LT. -180. ODO)  LONG  =  LONG  ♦  360. ODO 
IF  (LONG. 6T. 180. ODO)  LONG  >  LONG  -  360. ODO 


20 

21 

22 

23 

24 

23 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 


L2  ■  LONG  «  Y 

IF  (ZONE. LE. O.ODO)  60  TO  1  8Z0NE  NOT  YET  COMPUTED 

C 

C  CHECK  FOR  ZONE  OVERLAP  ON  EITHER  SIDE  OF  180/-180  DEGREES 
C  L0N6ITUDE.  FORCE  CONSISTENT  LONGITUDE  VALUE. 

C 

IF  (LONG. LT. O.ODO. AND. ZONE. EQ. 60. ODO)  LONG  =  LONG  +  360. ODO 
IF  (LONG. GT. O.ODO. AND. ZONE. EQ.1. ODO)  LONG  >  LONG  -  360. ODO 
L2  ■  LONG  *  Y  •RECOMPUTE  L2  FOR  OVERLAP 

C  E  IS  THE  ECCENTRICITY  SQUARED 
1  E*(A1-B1 ) * ( A 1 ♦ B 1 ) / A 1 / A 1 

C  A,  B,  C,  D  DEPEND  ONLY  ON  E  AND  ARE  AS  IN  THE  URITE  UP 
A*1  .♦£*( . 75+E*(0. 7031 25D0+E*. 68359375D0) ) 

B*E* (0.373+E*(0.46875+E»0. 51 29531 25D0>> 

C  *E • E • ( .05839375D0+E* . 1 025390625D0 ) 

DzE*E*E*O.Ot 1 3932291 6666667D0 
CI-DCOS(LI) 

S1*DSIN(L1) 

T«DTAN(L1 ) 

S2»C1*C1*E/(1.-E) 

T2»T*T 


A1  6T  B1 . 
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FLHUTM  2  of  2 


41 

C 

N  IS  THE  RADIUS  OF  CURVATURE  IN  THE  PRIME  VERTICAL 

42 

N*A1/DS0RT ( 1 .-E*S1 *S1 ) 

43 

C3«1 .-T2+S2 

44 

C4*T*(5.-T2+S2*(9.+4.*S2)> 

45 

C5*5.+T2*(T2-18. )*S2*(14.-58.*T2) 

44 

C4=T*(41 .  +T2*(T2-58. )+S2* (270.-330 ,*T2) ) 

47 

IF  (20NE.6T.0.0D0)  60  TO  3  IFORCE  PREVIOUS  ZONE 

46 

IF  (L2.6E.0.0)  60  TO  2 

49 

C 

2  IS  THE  ZONE 

50 

Z*30. ♦DINK  LONG/4. ) 

51 

IF  (Z.EQ.O.ODO)  Z  >  40.0D0  6F0R  -180  DEGREES  LONGITUDE 

52 

60  TO  3 

53 

2 

Z»31 ,+DINT (LONG/4. ) 

54 

C 

H  IS  THE  CENTRAL  MERIDIAN  IN  RADIANS 

55 

3 

N*(4.*Z-183. )*T 

54 

C 

S  IS  THE  ARC  LENGTH  ALONG  THE  CENTRAL  MERIDIAN 

57 

S«A1*(1.-E)*(A*L1-B*DSIN(2.*L1)+  C*DSIN(4.*L 1 )-D*DSIN(4,*L1 ) ) 

56 

D2=<L2-M)*C1 

59 

03>02*02 

40 

c 

CHECK  TO  AVOID  UNDERFLOU 

41 

IF  <A8S(D2) .LT. 1 . 0E-5 )  60  TO  10 

42 

E2s  N*D2*( 1 .♦03*<C3/4.+D3*C5/120. > )*.9994D0 

43 

N2=.9994D0»S+.9994D0*N*D3*(.5*T+D3*(C4/24.+D3*C4/72O.>) 

44 

60  TO  20 

45 

10 

E2-  N*D2*< 1 ,+D3*C3/4. >*.999400 

44 

N2*.9994D0*S+.9994P0*N*D3*T*.5 

47 

20 

E2*E2+5.D5 

48 

IF ( N2.6E .0. )  60  TO  30 

49 

N2»M2*t.D7 

70 

C 

0(1)  (OR  E2)  IS  THE  UTM  EASTING 

71 

C 

0(2)  (OR  M2 )  IS  THE  UTM  NORTHING 

72 

30 

0(  1  )»£2 

73 

0(2)*N2 

74 

IF  (Z0NE.6T.0.0D0)  60  TO  899 

75 

ZONE  «  2 

74 

899 

RETURN 

77 

END 
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10 

11 

12 

13 

14 

15 

16 

17 

18 
IV 
20 
21 
22 


SUBROUTINE  FLHXYZ(FI,X,AXES> 

C 

C  THIS  PROGRAM  TRANSFORMS  GEOGRAPHIC  COORDINATES  TO  LOCAL 
C 

IMPLICIT  DOUBLE  PRECISION  <A-H,0-Z> 

DIMENSION  FI(3),X(3),AXES(2> 

C 

SINFI  =  DSIN  (fl(D) 

cosfi  =  dcos  (Finn 

SINL  *  DSIN  <FI(2) ) 

COSL  =  DCOS  ( F I < 2 ) ) 

ALFA  *  (AXES(2)/AXES( 1 ) )**2 

E2  =  1.  -  ALFA 

6AHA  ■  1.  -  E2*(SINFI**2) 

6AMA  =  DSQRT  (GAMA) 

FN  =  AXES ( 1 ) /GAMA 

6AHA  =  (FN+FI (3) )  *  COSFI 

X ( 1 )  *  GAMA  *  COSL 

X ( 2 )  =  GAMA  *  SINL 

X ( 3 )  =  (FN  *  ALFA  +  FK3) )  ♦  SINFI 

RETURN 

END 


FLH00100 

FLH00200 

FLH00300 

FLH00400 

FLH00500 

FLH00600 

FLH00700 

FLHOO 

FLHOO 

FLHOO 

FLHOO 

FLH01200 

FLH01300 

FLH01400 

FLHO 

FLH01600 
FLH01700 
FLHO 1 800 
FLH01900 
FLH02000 
FLH021 00 
FLH02200 


1 

SUBROUTINE  LCGO (PIN , POUT, AXES, ORIG, HO, NODE ) 

LCG00100 

2 

C 

LCG00200 

3 

C 

THIS  PROGRAM  TRANSFORMS  COORDINATES  FROM  LOCAL  TO  GEOCENTRIC 

AND  VICE  LCG00300 

4 

C 

VERSA 

LCG00400 

3 

C 

LCG00500 

6 

C 

MODE  =  1  ,  GEOCENTRIC  TO  LOCAL 

LCG00600 

7 

C 

MODE  =  2  ,  LOCAL  TO  GEOCENTRIC 

LCG00700 

8 

C 

LCG00800 

V 

IMPLICIT  DOUBLE  PRECISION  <A-H,0-Z> 

LCGOOVOO 

10 

DIMENSION  FIN( 3) , POUT <3),AXES(2),0RIG(3),RT(3,3),P0(3> 

LCG01000 

ft 

C 

LCG01100 

12 

POUT  < 1 )  =  OR  I G  < 1 ) 

LCG01200 

13 

POUT ( 2 )  =  OR  I G ( 2 ) 

LCG01300 

14 

P0UH3)  =  HO 

LCG01400 

15 

CALL  FLHXYZ(POUT,PQ,AXES> 

LCGO 1500 

16 

CALL  0RTM6A(0RIS,RT) 

LCG01 600 

17 

60  TO  (10, 30), NODE 

LCG01700 

18 

10  DO  20  I  =  1,3 

LCG01800 

IV 

POUT ( I )  *  0.0 

LCGOIVOO 

20 

DO  20  J  *  1,3 

LCG02000 

21 

POUT ( I )  *  POUT ( I )  +RT ( I , J) *<PIN( J) -P0( J) ) 

LCG02100 

22 

20  CONTINUE 

LCG02200 

23 

60  TO  50 

LCG02300 

24 

30  DO  40  I  =  1,3 

LCG02400 

23 

POUT ( I )  *  POt  I ) 

LCG02500 

26 

DO  40  J  *  1,3 

LCG02600 

27 

POUT ( I )  *  POUT (I)*RT(J,I)*PIN(J) 

LCG02700 

28 

40  CONTINUE 

LC002800 

2V 

30  RETURN 

LCG02V00 

30 

END 

LCG03000 

A-  I  r> 


SUBROUTINE  MOVE<A,B,N) 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 


C 

C 

IMPLICIT  DOUBLE  PRECISION  <A-H,0-Z> 
C 
C 
C 
C 

DIMENSION  A<1),B<1) 

DO  10  I  =  1,N 
10  B1I)  =  Nil) 

RETURN 

END 


MOVOOIOO 
MOV00200 
NOV00300 
NQV00400 
MOV00500 
MOV00600 
HOV00700 
HOV00800 
HOV00900 
MOVOtOOO 
MO V0 1100 
NOVO 1200 
H0V01300 


1 

2 

SUBROUTINE  0RTM4A<V,A> 

ORTOOIOO 

ORT00200 

3 

IMPLICIT  DOUBLE  PRECISION  <A-H,0-Z) 

ORT00300 

4 

DIMENSION  V( 3) , A ( 3 , 3) , C 1 3) , S ( 3 ) 

0RT00400 

5 

C 

ORT00500 

6 

DO  10  I  *  1,3 

OR T 00600 

7 

C(I)  «  COS  (V(I1) 

ORT00700 

e 

10  StI)  =  SIN  (V(D) 

ORTOOSOO 

9 

A(1,1)  «  S(3)*S<1)*C<21  -  S<2>»C<3> 

0RT00900 

10 

A( 1 , 2)  *  S(3)»S(1 )*S(2)  ♦  C(2)*C(3) 

ORTOIOOO 

ii 

A( 1 ,3)  *  -S(3)*C<  1 ) 

0RT01 100 

12 

A(2,1)  *  -C(3)*S( 1 )*C<2)-S(2)*S(3) 

0RT0120C 

13 

A C 2 , 2 )  »  -C(3)*S<1)«S(2>  *  C(2)*S(3> 

ORT01300 

14 

A<2,3)  «  C(3)*C(1 1 

0RT01400 

15 

A(3, 1 )  «  C<1)*C<2) 

ORT01500 

n 

A(3,2)  =  C(1)*S<2) 

ORT01400 

17 

A(3,3)  *  S(1) 

OR TO 1700 

18 

RETURN 

ORT01800 

19 

END 

ORT01900 
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SUBROUTINE  TRANCD  <R,S, SEMI, P, HO, FLAG, ZONE, KFLA6) 


IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

DIMENSION  R<3),S(3),P(3),F(4) ,SEMI(2) ,RP<3> 

C 

C  KFLA6  DETERMINES  THE  FORM  OF  THE  OUTPUT  AND  INDICATES  THE  FORM  OF 
C  THE  INPUT.  THE  POSSIBILITIES  ARE  AS  FOLLOWS 
C  KFLA6-1  MEANS  UTM  TO  6E0GRAPHIC 
C  <2  MEANS  UTM  TO  6E0CENTRIC 

C  *3  MEANS  UTM  TO  LOCAL  RECTAN6ULAR 

C  >4  MEANS  UPS  TO  GEOGRAPHIC 

C  *5  MEANS  UPS  TO  6E0CENTRIC 

C  =4  MEANS  UPS  TO  LOCAL  RECTANGULAR 

C  »7  MEANS  GEOGRAPHIC  TO  UTM 

C  =8  MEANS  6E06RAPHIC  TO  UPS 

C  *9  MEANS  6E06RAPH1C  TO  6E0CENTRIC 

C  <10  MEANS  GEOGRAPHIC  TO  LOCAL  RECTANGULAR 

C  <11  MEANS  GEOCENTRIC  TO  UTM 

C  <12  MEANS  6E0CENTRIC  TO  UPS 

C  <13  MEANS  GEOCENTRIC  TO  6E06RAPHIC 

C  <14  MEANS  6E0CENTRIC  TO  LOCAL  RECTANGULAR 

C  <13  MEANS  LOCAL  RECTANGULAR  TO  GEOCENTRIC 

C  <14  MEANS  LOCAL  RECTAN6ULAR  TO  6E0GRAPHIC 


TRANCD 


part  1  of  2 


TRA001 00 
TRA00200 
TRA00300 
TRA00400 
TRA00500 
TRAQ0400 
TRA00700 
TRA00800 
TRA00900 
TRA01000 
TRA01100 
TRA01200 
TRA01300 
TRA01400 
TRA01500 
TRA01400 
TRA01700 
TRAOIBOO 
TRA01900 
TRA02000 
TRA02100 
TRA02200 
TRA02300 


C  >17  MEANS  LOCAL  RECTANGULAR  TO  UTM 

C  <18  MEANS  LOCAL  RECTANGULAR  TO  UPS 

C 

C  P(3X1)  IS  A  VECTOR,  (LATITUDE, LONGITUDE, AZIMUTH)  OF  THE  LOCAL  SYSTEM. 
C  THE  ALTITUDE  OF  THE  LOCAL  RECTANGULAR  SYSTEM  IS  ZERO. 

C  A  AND  B  ARE  THE  MAJOR  AND  MINOR  SEMIAXES  OF  THE  SPHEROID. 

C  FLAG  <  0.  IF  IN  THE  NORTHERN  HEMISPHERE  ,  <1  IN  THE  SOUTHERN 
C  R  AND  S  ARE  THE  INPUT  AND  OUTPUT  COORDINATES.  RESPECTIVELY 
C  E  X  LATITUDE 

C  S=  N  OR  «  Y  OR  *  LONGITUDE  SIMILARLY  FOR  R. 

C  H  Z  ALTITUDE 

C 


IF  (KFLA6)  220,220,10 
10  IF  (KFLA6-18)  20,20,220 
20  IF  (KFLA6-4)  100,100,25 
23  FLA6  <  0.0 

IF  (KFLAG-7)  30,40,30 
30  IF  (KFLA6-8)  40,40,40 
40  IF  (R< 1  ) )  90,100,100 
40  IF  (KFLAG-1 1 )  70,80,70 
70  IF  (KFLAS-12)  100,80,100 
80  IF  <R<3>)  90,100,100 
90  FLA6  <  1.0 
100  CALL  MOVE  <R,RP,3> 

GO  T0( 1 ,1 , 1 ,2,2, 2,3, 4, 3 ,3,4, 4, 4,7,8, 8,8,8) ,KFLA6 
C 

C  UTM  TO  GEOGRAPHIC 
C 


1  F(1)  •  «P(1) 

F(2)»RP(2> 

F(3)*Z0NE 
F(4)*8£HI ( 1 ) 

F (5><SEMI (2) 

FUXFLAG 

CALL  UTMFLH  (F,B> 

8(3)*RP(3) 

IF  (KFLAG-1 >  110,220,110 
110  CALL  M0VE(8,RP,3) 


TRA02400 

TRA02500 

TRA02400 

TRA02700 

TRA02800 

TRA02900 

TRA03000 

TRA03100 

TRA03200 

TRA03300 

TRA03400 

TRA03500 

TRA03400 

TRA03700 

TRA03800 

TRA03900 

TRA04000 

TRA04100 

TRA04200 

TRA04300 

TRA04400 

TRA04500 

TRA04400 

TRA04700 

TRA04800 

TRA04900 

TRA03000 

TRA03I00 

TRA05200 

TRA05300 

TRA05400 

TNA03500 

TRA05400 

TIA03700 

TRA03800 

TRA03900 

TRA04000 

TRA04100 
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43 

44 

43 

44 

47 

48 

49 

70 

71 

72 

73 

74 

75 
74 

77 

78 

79 

80 
81 
82 

83 

84 

85 
84 

87 

88 

89 

90 

91 

92 

93 

04 

f: 

94 

97 

98 

99 
100 
101 
102 

1 03 

104 

103 

104 

107 

108 

109 

110 
111 
112 

113 

114 

113 

114 

117 

118 

119 

120 
121 
122 

123 

124 

125 


C 

C  6EQ6RAPHIC  TO  GEOCENTRIC 
C 

3  COLL  FLMXTZ  <RP,S,SEHI> 

IF  (KFLA6-2)  120,220,120 
120  IF  (KFLA6-3)  130,220,130 
130  IF  (KFLA6-9)  140,220,140 
140  CALL  HOVE  <S,RP,3> 

C 

C  6E0CENTRIC  TO  LOCAL  RECTANGULAR 
C 

7  CALL  LC60  <RP,S,SENI,P,H0,1 ) 
80  TO  220 

C 

C  UPS  TO  6E06RAPH1C 
C 

2  F(1)«RP<1> 

F(2)«RP(2) 

F ( 3 )  *  FLAB 
F ( 4 )  =  SEMI ( 1 ) 

F(5)  *  SENK2) 

CALL  UPSFLH  (F,S) 

S<3)*RP(3) 

IF  (KFLA6-4)  110,220,110 
C 

C  GEOGRAPHIC  TO  UTH 
C 

3  F ( 1 )=RP( 1 ) 

F (2)'RP(2) 

F ( 3 )  «  SEMI ( 1 ) 

F ( 4 )  •  SEMI (2) 

CALL  FLMUTH  (F,S,ZONE) 
S(3)'RP(3) 

60  TO  220 
C 

C  GEOGRAPHIC  TO  UPS 
C 

4  F ( 1 )=RP< 1 ) 

F(2)=RP(2) 

F  <  3 )  «  FLAG 
F ( 4 )  *  SEMI < 1 ) 

F (5)  =  SEMI (2 ) 

CALL  FLHUPS  <F,S> 

S(3)*RP(3) 

80  TO  220 
C 

C  GEOCENTRIC  TO  GEOGRAPHIC 
C 

4  CALL  XTZFLH(RP,S,SEHI ) 

IF  (NFLA6-13)  170,220,170 
170  IF  (KFLAG-14)  180,220,180 
180  CALL  HOVE  <S,RP,3> 

IF  (NFLAG-11)  190,3,190 
190  IF  (KFLA6-12)  200,4,200 
200  IF  (KFLAG-17)  4,3.4 
C 

C  LOCAL  RECTANGULAR  TO  GEOCENTRIC 
C 

8  CALL  LC80  <RP,S,SEHI ,P,N0,2) 
CALL  HOVE  (S,RP,3) 

IF  (KFLAG-15)  4,220,4 
220  RETURN 
END 


TRANCD 
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TRA04200 
TRA04300 
TRA04400 
TRA04300 
TRA04400 
TRA04700 
TRA04800 
TRA04900 
TRA07000 
TRA07100 
TRA07200 
TRA07300 
TRA07400 
TRA07500 
TRA07400 
TRA07700 
TRA07800 
TRA07900 
TRA08000 
TRA08100 
TRA08200 
TRA08300 
TRA08400 
TRA085C0 
TRA08400 
TRA08700 
TRA08800 
TRA08900 
TRA09000 
TRA09100 
TRA09200 
TRA09300 
TRA09400 
TRA09500 
TRA09400 
TRA09700 
TRA09800 
TRA09900 
TRA10000 
TRA10100 
TRA10200 
TRA10300 
TRA10400 
TRA1Q500 
TRA10400 
TRA10700 
TRA10800 
TRA10900 
TRA11000 
TRA1 1 100 
TRA11200 
TRA1I300 
TRA11400 
TRA11300 
TRA11400 
TRA11700 
TRA11BOO 
TRA11900 
TRA12000 
TRA12100 
TRA12200 
TRA12300 
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SUBROUTINE  UPSFLH < PT IN f PTOUT ) 

C 

C  THIS  SUBROUTINE  TRANSFORMS  UNIVERSAL  POLAR  STEREOGRAPHIC  COORDINATES 
C  TO  GEOGRAPHIC. 

C 

IMPLICIT  DOUBLE  PRECISION  <A-H,Q-Z> 

DIMENSION  PTIN( 5) ,PTOUT ( 2 ) 

C 

PI  =  3.141592653589793 
CON  =  1. 

IF(PTIN(3) )  20,10,20 
10  T  -  2000000.  -  PTIN12) 

FACT  =1.0 
60  TO  30 

20  T  =  PTIN12)  -  2000000. 

FACT  =  -1.0 

30  X  =  PTIN(I)  -  2000000. 

PTOUT (2)  *  AT AN2 ( Y , X > 

IF  (ABS  (X)  -  ABS  (Y>)  40,40,50 
40  R  =  Y  /  SIN  (PTOUT (2) ) 

GO  TO  60 

50  R  «  X  /  COS  (PTOUT (2) ) 

60  A  =  PT IN ( 4  > 

B  =  PT IN  <  5 ) 

E2  =  1.  -<B/A>**2 
E  =  SORT  (E2) 

TEMP  *  (1.  -  E)  /  (1.  ♦  E) 

EX  =  E/2. 

TAN2  =  (R  /  A)  •  (B/A)  /((TEMPmEX)  *  CON) 

P  »  ATAN  ( TAN2 )  *  2.0 
DO  70  I  =  1,10 

51  =  SIN  (P/2.0) 

Cl  =  COS  (P/2.0) 

52  *  SIN  (P) 

C2  =  COS  (P) 

TEMP  =  (1.  ♦  E  *  C2)  /  (1.  -  E*C2) 

TANP  =  SI  /  Cl 
C  =  TANP  *(TEHP«*EX)  -  TANZ 
C12  «  C1M2 

CONST  =  (1.  -  E  *  C2>  **2 
TEMPI  =  TEMP  MEX 
EXX  =  EX  -  1. 

TEHP2  «  TEMP  **EXX 

CONST  «  (.5  ♦  TEHP1/C12)  -  <E2  *  TEHP2  »  TANP  *  S2  /  CONST) 

DP  «  -  C  /  CONST 
P  «  P  ♦  DP 

IF  (ABS  (BP)  -  .000000005)  80,80,70 
70  CONTINUE 
PRINT  1000 

1000  FORMAT  (14H2  ERROR  UPSFLH) 

80  PTOUT < 1 )  *  (PI/2.  -  P)  •  FACT 
RETURN 
END 


UPSOOIOO 

UPS00200 

UPS00300 

UPS00400 

UPS00500 

UPS00600 

UPS00700 

UPS00800 

UPS00900 

UPS01000 

UPS01100 

UPS01200 

UPS01300 

UPS01400 

UPS01500 

UPS01600 

UPS01700 

UPS01800 

UPS01900 

UPS02000 

UPS02100 

UPS02200 

UPS02300 

UPS02400 

UPS02500 

UPSC2600 

UPS02700 

UPS02800 

UPS02900 

UPS03000 

UPS03100 

UPS03200 

UPS03300 

UPS03400 

UPS03500 

UPS03600 

UPS03700 

UPS03800 

UPS03900 

UPS04000 

UPS04100 

UPS04200 

UPS04300 

UPS04400 

UPS04500 

UPS04600 

UPS04700 

UPS04800 

UPS04900 

UPS05000 

UPS05100 

UPS05200 

UPS05300 
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SUBROUTINE  UTNFLH(F,0) 

C  UTH  TO  6E06RAPHIC 

IMPLICIT  DOUBLE  PRECISION  (A-H.L-2) 

DIMENSION  F(4),0(2) 

C  F  AND  0  ARE  THE  SAME  CALLING  VARIABLES  THAT  APPEAR  IN  6ENC0R 
C  F< 1 )  (OR  El)  IS  THE  UTM  EASTING.  F(2>  (OR  N1>  IS  THE  UTH  N0RTHIN6 
C  F ( 3 )  (OR  Z)  IS  THE  ZONE 

C  F<4)  (OR  A1>  AND  F(3)  (OR  St)  ARE  THE  AXES  OF  THE  SPEROID.  A1  GT  B1. 

C  F ( 4 )  (OR  H)  IS  1  OR  0.  1  FOR  THE  SOUTHERN  HEMISPHERE.  0  FOR  THE 
C  NORTHERN  HEMISPHERE 
E1«F(1 ) 

N1«F(2> 

Z =F  <  3 ) 

AW14) 

B1=F(5) 

H=F(4> 

C  Y  IS  A  CONVERSION  FACTOR,  DEG  TO  RAD 
Y=. 01 74532925 1994 32958D0 
C  E  IS  THE  ECCENTRICITY  SQUARED 
E*( A1-B1 )*( A1 +81 )/A1/A1 

C  A,  B,  C,  D  DEPEND  ONLY  ON  E  AND  ARE  AS  IN  THE  URITE  UP 
Ar1 .+£*( . 75+E* (0.7031 25D0+£*.48359375D0) ) 

B=E*(0.375+E*(0.44875+E*0.51 29531 25D0  > ) 

C=E*E*( .05859375D0+E*. 1025390625DO) 

D«E*E*E*.01 1 3932291 444464447D0 
C  E2  AND  N2  ARE  THE  UNSCALED  UTHS 
E2*(E1 -5.D5) /0.9996DO 
N2«(N1-H*1.D7)/0.9994D0 
C  H  IS  THE  CENTRAL  MERIDIAN  IN  RADIANS 
H*(4.*Z-183. )*Y 

C  THE  NEXT  4  LINES  CONTAIN  THE  NEWTON  ITERATION  TO  COMPUTE  THE  FOOTPOINT  LAT 
0*N2/A1 /( 1 ,-E )/A 
DO  10  1=1,5 

6*N2-A1 *<  1  .-E )*( A*Q-B*DSIN(2.*Q)+C*DSIN( 4 . *Q)-4.*D’*DSIN(4.*Q ) ) 

G1*-A1»(1 .-E)*(A-2.«B*DC0S(2.*O>+4.*C*DC0S<4.*Q>-4.*D*DC0S(4.*0>> 
0*0-6/61 
10  CONTINUE 

C  THE  FINAL  VALUE  OF  Q  IS  THE  FOOTPOINT  LATITUDE 
C1=DC0S(0> 

S1«BSIN(6) 

T*DTAN(8) 

T2*T*T 

C  N  IS  THE  RADIUS  OF  CURVATURE  IN  THE  PRIME  VERTICAL  OF  THE  FOOTPOINT 
N*A1/DSQRT ( 1 ,-E*S1 »S1 ) 

S2-C1  *C1 *E/ ( 1  - -E ) 

D2*E2/N 

B3«D2»D2 

L2EB2*(1 . *D3*<-( 1 . ♦2.*T2+S2)/4.*D3+(5.+4.*S2+T2*(28.+8.*S2+24.*T2) 
2>/120. ) )/C1*N 

L1*8+D3«T*(1.+S2)«(-.5+D3*(5.+S2*(1.-4.*S2)+T2*(3.-9.*S2))/24.) 

C  0(1)  (OR  LI)  IS  THE  LATITUDE  IN  RAD 

C  0(2)  (OR  L2)  IS  THE  L0N6ITUDE  IN  RAD 

0( 1 )*L1 
LONS  *  L2/Y 

IF  (L0N6.LT. -180. ODO)  LONG  *  L0N6  ♦  340. ODO 
IF  (L0N6.6T. 180.060)  LONG  >  LONG  -  340. ODO 
L2  «  LONG  •  Y 
0(2>-L2 
RETURN 
END 


A-.’n 


L 

%  I  * _ i  _„1 


1 

SU8R0UTINE  XYZFLH(XYZ,FLN, 

AXES) 

XYZOOIOO 

2 

C 

XY200200 

3 

C 

THIS  SUBROUTINE  CONVERTS  POINT 

COORDINATES  FROM  GEOCENTRIC  TO 

XY200300 

4 

C 

6E06RAPHIC 

XYZ00400 

3 

C 

THE  SUBROUTINE  USES  DOUBLE  PRECISION  ARITHMATICS 

XYZ00500 

4 

C 

XYZ00400 

7 

c 

SPECIFICATIONS 

XYZ00700 

8 

c 

XYZ00800 

9 

IMPLICIT  DOUBLE  PRECISION 

<A-H,0-Z) 

XYZ00900 

10 

DIMENSION  XYZ(3),FLH(3),AXES(2) 

XYZ01000 

11 

PI  *3.14159245358979300 

XYZ010 

12 

c 

XYZ01200 

13 

X  =  XY2 ( 1 ) 

XY201300 

14 

Y  =  XYZ<2> 

XYZ01400 

15 

Z  *  X YZ  <  3 ) 

XYZ01500 

14 

c 

XYZ01400 

1? 

c 

COMPUTE  L0N6ITUDE  LAMDA 

XYZ01700 

18 

c 

XYZ01800 

19 

IF  (X)  50,10,40 

XYZ01900 

20 

10  IF  < Y )  20,30,40 

XYZ02000 

21 

20  FLH(2)  =  -PI/2.0 

XYZ021 00 

22 

80  TO  80 

XYZ02200 

23 

30  FLH<2)  *  0.0 

XYZ02300 

24 

IF  (Z)  32,32,33 

XYZ02400 

25 

32  FLH(1)  *  -PI/2.0 

XYZ02500 

24 

GO  TO  34 

XYZ02400 

27 

33  FLH(1)  =  PI/2.0 

XYZ02700 

28 

34  FLH(3)  =  A  -  ABS  (Z) 

XYZ02800 

29 

GO  TO  120 

XYZ02900 

30 

40  FLH < 2 )  *  PI/2.0 

XYZ03000 

31 

GO  TO  80 

XYZ03100 

32 

50  CON  «  -PI 

XYZ03200 

33 

GO  TO  70 

XYZ03300 

34 

40  CON  *  0.0 

XYZ03400 

35 

70  FLH < 2 )  =  DATAN  (Y/X)  ♦  CON 

XYZOO 

34 

c 

XYZ03600 

37 

c 

COMPUTE  LATITUDE  PHI 

XYZ03700 

38 

c 

XYZ03800 

39 

80  A  *  AXES < 1 ) 

XYZ03900 

40 

B  «  AXES < 2 ) 

XYZ04000 

41 

E2  *  1.  -  (B/A)««2 

XYZ04100 

42 

T1  «  E2  *  Z 

XYZ04200 

43 

DO  100  I  *  1,10 

XYZ04300 

44 

ZP  «  T1  ♦  Z 

XYZ04400 

43 

R  *  <X**2  ♦  Y«*2  ♦  ZP««2) 

XYZ04500 

44 

SINPI  ■  ZP  /  DSORT  (R) 

XYZOO 

47 

T2  *  (A  •  E2  *  SINPI)  /  DSORT  (1.0  -  E2  *  S!NPI**2> 

XYZOO 

48 

D  *  DABS  < T 1  -  T2) 

XYZOO 

49 

IF  (D-.005)  110,110,90 

XYZ04900 

50 

90  T1  «  T2 

XYZ05000 

51 

100  CONTINUE 

XYZ05100 

52 

WRITE  (4,1000) 

XYZ05200 

53 

1000  FORMAT  (14H2  ERROR  XYZFLH) 

XYZ05300 

54 

110  R  *  X*«2  ♦  Y»»2 

XYZ05400 

55 

R  -DSORT  (R) 

XYZ050 

54 

ZP  *  Z  ♦  T2 

XYZ05400 

57 

FLH ( 1 )  *  BATAN  (ZP/R) 

XYZOO 

38 

DN  *  A  /DSORT  (1.  -  E2  *  SINPI **2 ) 

XYZ050 

39 

R  *  ZP  /  SINPI 

XYZ05900 

40 

FLH(3)  •  R  -  BN 

XY204000 

41 

120  RETURN 

XYZ04100 

42 

END 

XYZ04200 

,\-2  I 
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QUOTATION 


QUOTATION  NO 

181830 

PHASE  REFER  TO  THIS  QUOTATION  NO 
IN  AIL  CORRESPONDENCE  AND  ORDERS 


.  LNK  corp. 

TO:  •  4321  Hartwick  Rd. 

•  College  Park,  MD  20740 

•  Attn:  Russ  Smith/Dr.  Kanal 


NEAREST  DGC  SALES  OFFICE 

7927  Jones  Branch  Dr. 
McLean,  VA  22102 

703/827-9600 


THANK  YOU  FOR  YOUR  INQUIRY  WE  ARE  PLEASED  TO  QUOTE  AS  FOLLOWS: 
DATE  REFERENCE  FREIGHT  CHARGES  TERMS  -! 

7/29/80  PREPAY  AND  ADO  ol*D.G  C 


FOB  POINT  OF  ORIGIN 


ITEM 

QUANTITY 

DESCRIPTION 

MTLY 

UNIT 

PRICE 

MAINT 

TOTAL 

PRICE 

UNIT 

LIST  PRICE 

DISC 

% 

UNIT 

NET  DRICE 

TOTAL 

DELUXE  S250 

(One  disk,  two  tares,  consol 

5) 

l 

1 

8635-NB 

S250  with  256KB 
memory  ( two-way 
interleaved) 

40,500 

15 

34,425.00 

34,425.00 

2 

1 

8641 

High  Speed  Hard¬ 
ware  floating  Pt. 

6,195 

15 

5,265.75 

5,265.75 

3 

1 

6026 

800/1600  Tape  w / 
Cont. 

15,500 

15 

13,175.00 

13,175.00 

4 

1 

6026 

800/1600  Add-on  Tap 
Drive 

o 

11,300 

15 

9,605.00 

9,605.00 

5 

1 

8650-B 

Dual  Peripheral  Bay 

» 

2,600 

NA 

2,600.00 

2,600.00 

6 

1 

6067-N 

50MB  Disk  with  Cont 

19,800 

15 

16,830.00 

16,830.00 

7 

1 

6040 

60  CPS  Terminal 
Printer 

2,650 

15 

2, 252.50 

2,252.50 

HARDWARE 

TOTAL - 

1 

84,153.25 

ATTACHMENTS: 

DGC  DikouM  AgfMmtnt  Form _ 

DGC  Progrom  bctnt*  Agr**m«nt  Form  500 
DGC  Progrom  Avoilobdity  Schedule  Form  501 
DGC  Mo«nt#no«K#  Cof>»rcK»  Form - 


Tt-i  »  f.  Srl  \J  9!  a*  '.i*.  r-Pp.*  *OC  3r  t'AVi  M-UM  TmE  LA?;  “'ri  :  • 

if-'  in  so  BY  !  ‘  *  A  v'vt^AL  CORPORATION  DGC  PRIOR  ?  C,  C  . r  •  :  IK- 

Wi-7  C-F*i-  THIS  C-w*A'.C‘N  AND  AN/  CPOPP  PLACED  AS  A  P-'.M  l-Ul  J’F  ! 
!  E  •  •.  ivjlVr .  V  TO  T*«t  TERMS  AND  CON  CHONS  HEREON  AND  CM  7*iE  : 

..r.  h  jl.-.-M  C-  Tkl  A .  T  ACHED  AC'K^INTS  ANY  CONTRAC  i  RESULTING  ?- 
q  .  iAti^N  m-jS7  BC  StONID  IN  SOUImBCRD.  MASS.  BY  A  DULY  AUTHOR  Til 
A’  VC  OF  LC-v  ANY  ADDiHONS,  MOC  HJATlCNS  OR  WAIVERS  OF  ANY  Of  T*--f  *•« 
Iv-NDiTION*.  CONTAINED  HEREIN  OR  ON  THE  ATTACHED  AGREEMENTS  SHALL  On  -  L 
1 1 v*  IF  IN  WRITING  AND  AGREED  TO  BY  AN  AUTHORIZED  REPRESENTATIVE  C?  L 
MARES  NO  V.  APR  AN'.  Its,  EITHER  EXPRtSS  OR  IMPLIED,  AS  TO  THE  fUNESS  OF  THE  F- 
OUo’cO  F» -'•>  ANY  PAP.UCUIAR  PURPOSE,  SUCH  DETERMINATION  TO  BE  THE  SOLE  R 
6IHTY  OF  THE  CUSTOMER. 


Zuckerman,  Federal  Acct.  Mgr. 

CUSTOMER  COPY 


*/• 


DataGeneral 

I*  Rout*-  9.  9*  »ttboro .  MuwchuMlu  OlbBl 
Ttl.phon*  617  366  *911 


QUOTATION 
CONTINUATION  SHEET 


quotation  no 

181830 _ 

THE ASE  REFER  TO  THIS  QUOTATION  nO 
IN  All  CORRESPONDENCE  AND  ORDERS 


7/29/80 


ITEM  QUANTITY 


UNIT 

LIST  PRICE 


UNIT 

NET  PRICE 


0  ADD  BURST  MULTIPLEXOR  CHA 


6067-HN 

50  MB  Disk  with  Cont 
(substitute  this 
model  for  6067-N  on 
previous  page) 

8642 

Burst  Multiplexor 
Channel  (add  this 
model  to  previous 
page) 

REVISED 

TOTAL - 

21,800  15  18,530.00!  18,530.00 


3,150  15  2,677.50}  2,677.50 


■  88,530.75 


•Model/Feature  Defined  in  Dota  General  Price  List 


B-2 

CUSTOMER  COPY 


See  Page  I  for  Signature 


7/29/80 


unit 

list  MICE 


UPER  DELUXE 


0  ADD  INTEGRAL  ARRAY  PROCESS 


652B  I/O  only  module 

8652-C  I/O  system  Module 

8315  I/O  Bus  Repeater 

8644  Integral  Array  Proce 

8638  Writeable  Control  st 


SVISED  TOTAL- 


unit 

NET  MICE 


765.00 


15  765.00; 


765.00 


765.00 


15  1,020.00;  1,020.00 
15  12,405.75  12,405.75 
15  3,570.00  3,570.00 


I  107,056.50 


Model/ reoture  Defined  in  Dofo  Generol  Price  List 


B-3 

CUSTOMER  COPY 


See  PoQe  1  for  Signoture 


Model/Feature  Defined  in  Data  General  Price  List 


See  Page  1  for  Signature 
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A  U  THQRiZED  REPRESEN  TA  77  VE 


PRIZED  RE  PRESENT  A  TtVE 


PERKIN-ELMER 


QUOTATION 


DIGITAL  EQUIPMENT  CORPORATION 

PHONE  AC  617  897-5111  TWX  710-347-0212-CABLE  DIGITAL  MAYN,  TELEX  94-84-57 


datf  22  September-  1980 
REFERENCE  GS-OOC-Q1892 

DISCOUNT  AGREEMENT  NO  _ 


PRODUCT  LINE  CODE 


r 

TO 

L_ 


L  nk  Corporation 
1*321  Hartwick  Rd 
Suite  321 

College  Park,  KD  2071-0 
Attn:  Russ  Smith 


NEAREST  DIGITAL  SALES  OFFICE 

8301  Professional  Place 
handover,  Maryland  2C765 


J 


THANK  YOU  FOR  YOUR  INQUIRY,  WE  ARE  PLEASEO  TO  QUOTE  AS  FOLLOWS 


ITEM 

QTV 

MODEL  NO 

DESCRIPTION 

V 

M0*(!hi  * 
Mt'N’ 

Apr  ufc 

*Ni  ’ 

DiV 

* 

UNIT  PRICE 

NET  AMOUNT 

1 

1 

11X1*  L-CA 

256KB  ECC  MCE,  K96L2, 

fll*2 

$ 

27,700 

TUJ8,  120V 

2 

LA3S-HA 

DECwriter  II 

5  16 

$ 

1,700 

3. 

1 

FP11-F 

FPP  for  ll** a 

6  16 

$ 

3,100 

h . 

T_ 

P.JM02-AA 

67ME  Disk  &  Controller 

5170 

$ 

25,700 

5- 

1 

RM02-AA 

67MB  Disk  Drive 

51 1*0 

$ 

19,300 

6. 

i 

TJE16-EA 

800/1600  Dpi  1*5  ips  tape  dr. 

5ll*7 

$ 

20,200 

7. 

1 

TE16-AE 

800/1600  bpi,l*5ips  tape  drive 

510L 

$ 

12,800 

6. 

1 

H9642-DB 

Expansion  Cabinet  without 

1/C 

$ 

1,300 

end  panels 

9 

1 

BA11-KV,' 

Expansion  Box 

1/C 

$ 

3,200 

10 

1 

BC03M-25 

Null  Modem  Cable  1 

5/C 

$ 

60 

11 

1 

DL11-WE 

EIA  Async  Line  RT  Clock 

5  6 

$ 

820 

12 

1 

QJ737-AD 

RSX11-M 

5210 

$ 

7,800 

13 

1 

QP230-AD 

Fortran  IV/RSX  11K 

$  35 

$ 

1,000 

ll* 

1 

DD11-DK 

DD11-D,  2-SU  for  BA11-K 

1/C 

$ 

860 

Subtotal : 

$125,51*0 

Less  GSA  Discount 

$  18,831 

$106,709 

SUB  ' 

hai.  fpma.n  t,rm  for  fiO  days  I**  datn  ufosr-  **n  -  *'  i  "C  t .  P  a  M  t  q..‘ 

(  Ci.p,..a 

,  r  p.  r* 

:: 

$106,709.00 

Stanflar 

tpr^-s 

and  condition*  0<  sale  appealing  on  '  e*r  er  se  heteo*  *nd  O'  inp  Ipte  i,  a*,  nolpfl  at*nir** 

•no 

StTJ*(**«*d 

PLUS 

t  •  So* 

**  I'P  rc 

nsuitmq  2  *  Computer  Speoa1  3*  DEC  *  •  L>-q  ta>  C •• 

•  P 

sc  rii"'  ao 

,M„P, . 

NST  ALLATlON 

O-  7 

RS-d  «*r! 

p"nt.  Systems  te 

t *•*  pi.-enase  and  <W  .r  »n<* 

1  ease- 

*t  *•* 

INSURANCE 

*»>  A  ■■ , 

COn'^aCF  ’nSu't-nq  IFlP 

Quniat*o'>  ti»*  at cn pted  *•  Uqi*1  a  :c-*p.  -ve  n",**  ?■»  .« 

.  A  •* 

•  fed  et  •esn-t  v  4« 

NET 

n*  0  1  ’ 

f  r 

«t'H  Pe  prowdert  o»"  p'ope»t»  «»»•  1*  .»•  t*ar-\  •  a-d  a  -  *  *  **C 

P*1- 

S  •  on  or  1 

eq  .  p"-f  1 

total  amount 

•  ~-3i1p  O’-IPSA  ■'5l'uCtiO',S 

to  the  contra1*  »’*■  0'J"i  Oa'nrt  ’  -  r  «.v  #•  -•  ••'f  P  ■■  •• 

II-  I  I 

EN  try- 


QUOTATION 


DIGITAL  EQUIPMENT  CORPORATION 

PHONE  AC  617  897-5111  TWX  710-347-0212-CABLE  DIGITAL  MAYN.  TELEX  94-84-57 


Quotation  numbep 

C-10Ul0U-806l _ 

— reCSala-lW  IH  I 

DATE  12  Sppt.Pmh^T-  1980 
REFERENCE _ 


DISCOUNT  AGREEMENT  NO 


PRODUCT  LINE  CODE 

r 

LHK  Corporation 
-321  Hartwick  Rd. 

Suite  321 

College  Park,  MD  207L0 
| _  Attn:  Russ  Smith 


NEAREST  DIGITAL  SALES  OFFICE 

n 

8301  Professional  Place 
Landover,  Maryland  20785 


J 


THANK  YOU  FOR  YOUR  INQUIRY.  WE  ARE  PLEASED  TO  QUOTE  AS  FOLLOWS 


■^1 

QTY 

MODEL  NO 

DESCRIPTION 

T,D‘ 

MOA’-i  1 

MAN- 

ADD  0* 

■V  •  -IMG 

5. 

UNIT  PRICE 

net  amount 

11 •  3  ^  A- Y  ? 

11/3-A  CPU,  256KB  MOS  Memory; 

$15,300 

$15,300 

Serial  Line  Interface 

z 

1 

rr  a- “A 

Floating  Point  Processor 

$  2,900 

$  2,900 

3 

z 

FvTl-Ot:— AA 

67ME  Disk  &  Controller 

$21*.  000 

$L8 ,000 

1 

Ll/OO— HA 

DECvriter  II  -  RS232 

$  1,700 

$  1,700 

5 

X 

TCE16-EA 

800/16C0  bpi  15  ips  tape  drive 

$20,200 

$20,200 

&  controller. 

6 

i 

TE16-AE 

800/160C  BPI  tape  Drive 

$12,800 

$12,800 

I 

k96o-ca 

Cabinet 

$  1,575 

$  1,575 

8 

PC03M-23 

Null  Modem  Cable 

$  60 

$  60 

9 

2 

BAli-KE 

Expansion  Box 

$  3,000 

$  6,000 

10 

1 

•V737-A1 

RSX11M 

$  7,800 

$  7,800 

11 

1 

...P230-A1 

FORTRAI!  IV/RSX11M 

$  1,000 

$  1,000 

TOTAL: 

$117,335.00 

(GSA)  155  Discount 

$  17,600. 25 

$  99,73^.75 

SUB  TOTAL 

’  ' 1  • 

-em,,.r  lr»«  K'  0.1s 

«f  «*  tt'f  (iUtt-  „r|pr,j,  KHort.lrOt*  'H  »*•  I.OJ  b*  D,q',',l  COu'Pmpn 

CO'POta 

"Z  pun' 

st.r  .).«• 

.vd  o'  san  ape 

-.f  **q  me  ••»***  se  I'fpi,*  arid  o'  tp'r»-A  as  noted  above 

and 

attached  hereto 

PLUS 

.  #  s  , 

6  *  Discount  aq'PPmenl 

INSTALLATION 

•A  • 

SvStpm,.  ••• 

»••••  pu'cease  and  Group  iphps  and 

between  purchase' 

PLUS 

’ 

aod  1  "d  1 

•'  above 

INSURANCE 

1-  •  * 

.  A  • 

••  *  .  ■  j  •• « 

a,,,  'v  •••  be  accepted  at  O-p-ta' s  co'OP'JtF  offices  bv  a  do 

autbon/ect  leb'eseoiat 

VP 

NET 

f  ;  .  1 

■  •:  . . . 

WI  be  C*' f  »v  deft  r.r  proper?*  V*h'le  in  bans-?  and  a  cha'qe  o’  $  60 

pe* 

S  tOP  00  O*  equipment 

TOTAL  AMOUNT 

«,•  ire  r.-i'?"!',  .jrf.  clPa-',  slated  or  tup  lace  0*  tre  P„'cnaSe'S 

O'dr 

ESTIMATED  DELIVERY  SCHEDULE  s-.'H  >i  »  to  mooiucation  p»  digital  , 


QUOTATION  PREPARED  BY 


Dale  Vogel 
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